{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vast 汽车价格离群值检测"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "离群点（Outliers），简单而言就是离其余数据点非常远的数据点。它们会极大的影响后续的分析结果，甚至产生有误导的分析结果。\n",
    "\n",
    "[Vast](http://www.vast.com/)向3个行业的出版商、市场和搜索引擎提供数据，这三个行业包括汽车、房地产和休闲、住宿和旅游。Vast的系统通过白标签集成，并在一些非常受欢迎的消费应用程序（Southwest GetAway Finder，AOL Travel，Yahoo! Travel，Car and Driver等等）中提供搜索结果、产品建议和特别优惠。\n",
    "\n",
    "Vast的汽车数据是由成千上万的二手汽车卖家提供，并且向市场公布。由于这些数据是由用户手工录入，因此容易受到人为失误的影响，比如用户在错误的字段中提交值，或者无意中发生错误或胖手指值。 对于8岁的车辆，里程表读数为100,000英里。 直觉告诉我们，100,000美元是大多数小型车的不寻常的价格。 一个上市的42000美元是合理的，比如2013年的凯迪拉克ATS豪华版，而对另一个（例如1997年的别克莱斯布雷）来说，这可能是意料之外的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![2013cadillac](./img/2013cadillac.png)\n",
    "![2013cadillac](./img/1997buick.png)\n",
    "<div align=\"center\">图1： 2013年的凯迪拉克ATS豪华版与1997年的别克莱斯布雷汽车信息</div>\n",
    "<br>\n",
    "检测离群值有利于纠正错误，向用户提供卓越和合适的产品。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. 数据集描述"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "数据集包括训练数据集和测试数据集：\n",
    "\n",
    "+ `accord_sedan_testing.csv`\n",
    "+ `accord_sedan_training.csv`\n",
    "\n",
    "![hondaaccord](./img/2016-honda-accord-sdn_100522223_l.jpg)\n",
    "测试数据集`accord_sedan_training.csv`包含了417辆本田雅阁轿车的信息列表。各个特征如下：\n",
    "\n",
    "|特征名称|含义|类型|取值示例|\n",
    "|:-----:|:------:|:----:|:----:|\n",
    "|$price$|价格|int|14995|\n",
    "|$mileage$|已行驶英里数|int|67697|\n",
    "|$year$|上市年份|int|2006|\n",
    "|$trim$|档次|str|ex;：高配款且带皮革内饰<br>ex：高配款<br>lx：低配款|\n",
    "|$engine$|引擎缸数|int|4 Cyl：4缸<br>6 Cyl：6缸|\n",
    "|$transmission$|换挡方式|str|Manual：手动挡<br>Automatic：自动挡|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. 导入数据集并切分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/envs/py27/lib/python2.7/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated since IPython 4.0. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.\n",
      "  \"`IPython.html.widgets` has moved to `ipywidgets`.\", ShimWarning)\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd # 读取数据表并进行基于DataFrame结构的操作\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn import preprocessing\n",
    "from sklearn.feature_extraction import DictVectorizer\n",
    "from sklearn.linear_model import LinearRegression\n",
    "%matplotlib inline\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore') # 不显示warning信息\n",
    "pd.options.display.width = 900 # Dataframe 显示宽度设置"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 画图配色方案\n",
    "+ 绿色 #5dbe80\n",
    "+ 蓝色 #2d9ed8\n",
    "+ 紫色 #a290c4\n",
    "+ 橙色 #efab40\n",
    "+ 红色 #EE5150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The shape is (417, 6)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>price</th>\n",
       "      <th>mileage</th>\n",
       "      <th>year</th>\n",
       "      <th>trim</th>\n",
       "      <th>engine</th>\n",
       "      <th>transmission</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>14995</td>\n",
       "      <td>67697</td>\n",
       "      <td>2006</td>\n",
       "      <td>ex</td>\n",
       "      <td>4 Cyl</td>\n",
       "      <td>Manual</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>11988</td>\n",
       "      <td>73738</td>\n",
       "      <td>2006</td>\n",
       "      <td>ex</td>\n",
       "      <td>4 Cyl</td>\n",
       "      <td>Manual</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>11999</td>\n",
       "      <td>80313</td>\n",
       "      <td>2006</td>\n",
       "      <td>lx</td>\n",
       "      <td>4 Cyl</td>\n",
       "      <td>Automatic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>12995</td>\n",
       "      <td>86096</td>\n",
       "      <td>2006</td>\n",
       "      <td>lx</td>\n",
       "      <td>4 Cyl</td>\n",
       "      <td>Automatic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>11333</td>\n",
       "      <td>79607</td>\n",
       "      <td>2006</td>\n",
       "      <td>lx</td>\n",
       "      <td>4 Cyl</td>\n",
       "      <td>Automatic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>10067</td>\n",
       "      <td>96966</td>\n",
       "      <td>2006</td>\n",
       "      <td>lx</td>\n",
       "      <td>4 Cyl</td>\n",
       "      <td>Automatic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>8999</td>\n",
       "      <td>126150</td>\n",
       "      <td>2006</td>\n",
       "      <td>lx</td>\n",
       "      <td>4 Cyl</td>\n",
       "      <td>Automatic</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   price  mileage  year trim engine transmission\n",
       "0  14995    67697  2006   ex  4 Cyl       Manual\n",
       "1  11988    73738  2006   ex  4 Cyl       Manual\n",
       "2  11999    80313  2006   lx  4 Cyl    Automatic\n",
       "3  12995    86096  2006   lx  4 Cyl    Automatic\n",
       "4  11333    79607  2006   lx  4 Cyl    Automatic\n",
       "5  10067    96966  2006   lx  4 Cyl    Automatic\n",
       "6   8999   126150  2006   lx  4 Cyl    Automatic"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train = pd.read_csv('accord_sedan_training.csv')\n",
    "print 'The shape is',train.shape\n",
    "train.head(7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 将 train 划分为 x 和 y\n",
    "train_x = train.drop('price',1)\n",
    "train_y = train['price']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. 特征提取并构建线性回归模型\n",
    "使用`sklearn.feature_extraction`中的`DictVectorizer`将名义型变量实现One-Hot编码，得到一个10维向量空间，包括：\n",
    "\n",
    "+ engine=4 Cyl\n",
    "+ engine=6 Cyl\n",
    "+ mileage\n",
    "+ price\n",
    "+ transmission=Automatic\n",
    "+ transmission=Manual\n",
    "+ trim=ex\n",
    "+ trim=exl\n",
    "+ trim=lx\n",
    "+ year\n",
    "\n",
    "注意：在异常值检测时需要利用到$~price~$特征。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimension is 10\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['engine=4 Cyl',\n",
       " 'engine=6 Cyl',\n",
       " 'mileage',\n",
       " 'price',\n",
       " 'transmission=Automatic',\n",
       " 'transmission=Manual',\n",
       " 'trim=ex',\n",
       " 'trim=exl',\n",
       " 'trim=lx',\n",
       " 'year']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# One-Hot编码方法一：使用`DictVectorizer`实现\n",
    "dv = DictVectorizer()\n",
    "dv.fit(train.T.to_dict().values()) # one-hot 编码\n",
    "print 'Dimension is',len(dv.feature_names_)\n",
    "dv.feature_names_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "也可以使用`pandas`的`get_dummies`函数实现one-hot编码："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mileage</th>\n",
       "      <th>year</th>\n",
       "      <th>engine_4 Cyl</th>\n",
       "      <th>engine_6 Cyl</th>\n",
       "      <th>trim_ex</th>\n",
       "      <th>trim_exl</th>\n",
       "      <th>trim_lx</th>\n",
       "      <th>transmission_Automatic</th>\n",
       "      <th>transmission_Manual</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>67697</td>\n",
       "      <td>2006</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>73738</td>\n",
       "      <td>2006</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>80313</td>\n",
       "      <td>2006</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>86096</td>\n",
       "      <td>2006</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>79607</td>\n",
       "      <td>2006</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   mileage  year  engine_4 Cyl  engine_6 Cyl  trim_ex  trim_exl  trim_lx  transmission_Automatic  transmission_Manual\n",
       "0    67697  2006             1             0        1         0        0                       0                    1\n",
       "1    73738  2006             1             0        1         0        0                       0                    1\n",
       "2    80313  2006             1             0        0         0        1                       1                    0\n",
       "3    86096  2006             1             0        0         0        1                       1                    0\n",
       "4    79607  2006             1             0        0         0        1                       1                    0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# One-Hot编码方法二：使用`pandas`的`get_dummies`函数实现\n",
    "nomial_var = ['engine','trim','transmission']\n",
    "multi_dummies = [] # 存储三个 DataFrame\n",
    "train_x_dummies = train_x[['mileage','year']]\n",
    "for col in nomial_var:\n",
    "    dummies = pd.get_dummies(train_x[col], prefix = col)\n",
    "    train_x_dummies = pd.concat([train_x_dummies,dummies],axis=1) # 将编码结果与非编码特征水平拼接起来\n",
    "train_x_dummies.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'12084.24 + (-337.20 engine=4 Cyl) + (337.20 engine=6 Cyl) + (-0.05 mileage) + (0.00 price) + (420.68 transmission=Automatic) + (-420.67 transmission=Manual) + (208.93 trim=ex) + (674.60 trim=exl) + (-883.53 trim=lx) + (2.23 year)'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 构建线性回归模型\n",
    "train_x_array = dv.transform(train_x.T.to_dict().values())\n",
    "# train_x_array = train_x_dummies.values # 也可以使用get_dummies得到的结果\n",
    "\n",
    "LR = LinearRegression().fit(train_x_array,train_y)\n",
    "' + '.join([format(LR.intercept_, '0.2f')] \n",
    "           + map(lambda (f,c): \"(%0.2f %s)\" % (c, f), \n",
    "               zip(dv.feature_names_, LR.coef_)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "即拟合得到的模型为：\n",
    "\n",
    "$price≈12084.24 − 337.20(engine=4Cyl) + 337.20(engine=6Cyl) - 0.05(mileage) + 420.68(transmission=Automatic) - 420.67(transmission=Manual) + 208.93(trim=ex) + 674.60(trim=exl) - 883.53(trim=lx) + 2.23(year)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. 离群值检测\n",
    "由测试集拟合得到的模型，我们可以预测测试集中的价格，计算每个样本的绝对误差，并得出"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1391.71504177,  2200.19252869,  2626.93955514,  3857.45997116])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred_y = LR.predict(dv.transform(train.T.to_dict().values()))\n",
    "train_error = abs(pred_y - train_y) # 计算绝对误差\n",
    "np.percentile(train_error,[75,90,95,99]) # 计算绝对误差数据的百分位数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们也可以画出 train_error 的盒图来观察离群值。\n",
    "![boxplot](./img/boxplot.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x107f5d950>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAEFCAYAAABAVTQtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADLNJREFUeJzt3X+s3Xddx/HXve3arV0pM9yZKFOIGx9DINRpnMqAxgCr\nAqshMzEGGRAxJPxhE5OpC8usAeMvEHUaCcomMP9xDBeXsC2KE/DHEEcDhOWzjaCu8McKtF1d13Zr\n6x/n3HLXnttd7+55n3O6xyNZ1vO9n34/38+5O8/z3ffe+71zJ0+eDAA15id9AADPJaILUEh0AQqJ\nLkAh0QUotP5sH9y379Cqv7Xhoos2Zf/+w6v961PDOqaLdUyXc2Ed41jDwsKWueU+NrYz3fXr141r\n16WsY7pYx3Q5F9ZRvQaXFwAKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoA\nhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSh01l9MOY0++MHfz8GDB9ZkX4cPP54k2bRp87Jj5ufn\ncuLEqn8/5xm2bn1+du26bs32B8yWmYvuwYMH8p3938m6TRuf9b6OHz2SJHmy6PfSHT98tGYiYGrN\nXHSTZN2mjbn4zduf9X4evf3eJFmTff1/5gOeu1zTBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKiS5A\nIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCoku\nQCHRBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJ\nLkAh0QUotH4cO73jjttywQUb8vrXXz2O3XMOuOOO25IkO3deM+EjgVpjOdPds+f+3HfffePYNeeI\nPXvuz54990/6MKCcywsAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSgkugCFRBegkOgCFBJdgEKi\nC1BIdAEKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBC\nogtQSHQBCokuQCHRBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKrZ/0AcAoDz3UT/15795Hsm/fo1lY\nuDgve1nLgQOHc9llbeT4lW4/fczevY+cevzCF17ytPGr3fdK5p5WS5//ZHJrmNRzOM55RZepdNdd\nd5768969j+To0aPZuHFjev9Knnzy+BkvhsXxK91++pizRXe1+17J3NNq6fOfTG4Nk3oOxzmv6DJ1\nHnqo5+GHHzxj+5EjT+SBBx44NWbxBbF0/Eq2P9NcDz/84Knxq933SuaeVqOek0msYVLP4bjnHUt0\nDx9+PMeOHcvu3dev+b4PHNifk+vm1ny/FU4cezIHntg/luflmczPz+XEiZPl8y7nwIH92bBhw8iP\nnX6WtdyYxRfD0vEr2b6SuRbHr3bfK5l7Wo16Tiaxhkk9h+Oe1xfSAAqN5Ux306bNufDCC3PDDe9b\n833v3n19Dh59fM33W2F+w3nZunFzbrzxd8rnXljYkn37DpXPu5yzne3v2PHG3HTTB87693fseOPI\n8SvZvpK5Fsevdt8rmXtajXpOJrGGST2H457XNV2mzmWXtVx66UtOPV76hbQXv/hFZ3whben4lWwf\nNddyX0hb7b5XMve0Ov35X9w2yeOonH/c84ouU2npGcaobxk72/iVbD99zOnRXYt9z9oZ7lLTcuyT\nOo5xziu6TKXTzyoXLXeZ5GxnsiuZ62zjVrvvWTvDXWpajn1SxzHOeX0hDaCQ6AIUEl2AQqILUEh0\nAQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSgkugCFRBegkOgCFBJdgEKiC1BI\ndAEKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQ\nSHQBCokuQCHRBSi0fhw73bbt8lxwwYZx7JpzxLZtl0/6EGAixhLdnTuvycLCluzbd2gcu+ccsHPn\nNZM+BJgIlxcACokuQCHRBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKiS5AIdEFKCS6AIVEF6CQ6AIU\nEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSgkugCFRBegkOgC\nFBJdgEKiC1BIdAEKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUGj9pA9gNY4fPppHb793DfZzJEnW\nZF8rm+9osnFzyVzAdJq56G7d+vw129fh44N/bzpLCOfn53LixMm1mXDj5jU9fmD2zFx0d+26rnS+\nhYUt2bfvUOmcwLnLNV2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHR\nBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKiS5AobmTJ9foN90C8Iyc6QIUEl2AQqILUEh0AQqJLkAh\n0QUoJLoAhdav9Q5ba/NJ/jzJK5IcTfLLvfeH13qetdBauyLJ7/Xet7fWLk1yS5KTSb6S5N299xOt\ntRuTvCHJU0l29d4/v9zYCRz/eUk+kuRFSTYmeW+Sr446tilfx7okH07ShsfyriRHZm0di1prFyf5\nzySvGx7nGcc27etord2f5LHhw68n+VCSP87geO/pve9e7rXeWvuJ08eWL2CotfabSa5OsmF4rP+c\nCX8+xnGm+3NJzu+9/2SS30jy/jHM8ay11q5L8pdJzh9u+kCS9/TeX5VkLsnO1trlSV6T5Iokv5Dk\nz5YbW3nsS7wlybeHx7EjyU2jjm0G1vGmJOm9vzLJe5K8b9SxzcA6Ft8IP5TkieWObdrX0Vo7P8lc\n73378J+3J/mLJL+Y5MokV7TWfiTLv9ZHjS3XWtue5KeSvDKD5/uSTMHnYxzRvTLJXUnSe//3JD82\nhjnWwteSvHnJ4x/N4F0wST6V5LUZrOWe3vvJ3vv/JFnfWltYZuwk/G2SG4Z/nsvgXXrm1tF7/7sk\nvzJ8+INJDmQG1zH0hxlE55vDx7O4jlck2dRau6e19unW2quTbOy9f633fjLJ3fnuOp72Wm+tPW+Z\nsZNwVZIvJ/lkkr9Pcmem4PMxjug+L8nBJY+Pt9bW/DLGs9V7/0SSJ5dsmhv+R5Ikh5JszZlrWdw+\namy53vv/9t4Ptda2JLktg7PEmVtHkvTen2qt/XWSP01y6zLHNtXraK29Lcm+3vvdSzbP3DqSHM7g\nzeOqDC713Dzctmi5dRwfbntsxNhJeEEGJ30/n8E6bk0yP+nPxzii+1iSLUvn6L0/NYZ51trSazVb\nMjjbOn0ti9tHjZ2I1tolSf4pycd673+TGV1HkvTer03ykgyu716w5EOzso53JHlda+3eJNuSfDTJ\nxUs+PivreDDJx4dnfg9mEKTvWfLx5dYxP2LbJNfx7SR3996P9d57Bl8nWBrOiXw+xhHdf0nys0ky\nvKD+5THMMQ5fHF4DSpKfSfLZDNZyVWttvrX2Axm8gXxrmbHlWmvfm+SeJL/ee//IcPMsruOXhl/w\nSAZnVCeSfGHW1tF7f3Xv/TW99+1J9iR5a5JPzdo6MnjzeH+StNa+L8mmJI+31n6otTaXwRnw4jqe\n9lrvvT+W5NiIsZPwuSQ7Wmtzw3VsTvKPk/58jON/+z+Zwbv9v2ZwnfHtY5hjHH4tyYdbaxuSPJDk\ntt778dbaZ5P8WwZvUO9ebuwkDjjJ9UkuSnJDa23x2u6vJvmTGVvH7Ulubq19Jsl5SXYNj2fWPh+j\nzOJ/V3+V5JbW2ucy+Mr9OzJ4I7w1yboMrn/e11r7j4x+rb/r9LHVC0iS3vudw+vRn893n+evZ8Kf\nD7d2BCjkhyMACokuQCHRBSgkugCFRBegkOgyU1prV7fWfnvSxwGr5VvGAApN3T0ReO4a/vTP7gzu\niXFJBt/U/t4kdyT5VgY/xvnxJNt7729rrb02g5+cmk/y3xnc2erxJH+QZHsG35x/S+/9j0oXAmfh\n8gLT5scz+ImgH87gtptvyOA+u2/pvZ+6y1NrbWMGP/V0be/95Um+lOTaJO9Mkt775cN97Wytvap0\nBXAWznSZNp8Z3pwkrbWPZXDLx0d77/912riXJ/lG731PkvTerx/+nduSbGut/fRw3IXDsZP6+X94\nGtFl2iy9I9388PETI8YtvS1nWmtbM7gT1Lok1/Xebx9uf0EGlxxgKri8wLS5srX2/cNfBfPWDG4e\nPUpPstBae+nw8XUZ3Gjl00ne2Vo7r7V2YQZ3mrpi3AcNKyW6TJtvZnAf2q8m+UaSfxg1qPd+JINf\nV/TR1tqXkrw0ye9m8FsbHkryxSRfSHJz7/3e8R82rIxvGWNqDL974beG96OFc5IzXYBCznQBCjnT\nBSgkugCFRBegkOgCFBJdgEL/BwwOy/d/V0n0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107f5dfd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.boxplot(x = train_error,palette = \"Set2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到测试集中存在部分离群值。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在本案例中，我们设定置信水平为0.95，即认为超过95%百分位数的train_error为离群值。下面我们在二维空间中画出正常值（蓝色）与离群值（红色）："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAckAAAGsCAYAAAC2BlyYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXucE9X5/z9J9sruArtsQLEqSvVoa9ufLeKtWtpvL19t\na1vsFy+12hu9WS9oa/vVVsF6qdpitbW2XqooX1CsKKiAF1BEFhTwAniZXRZ2gb2w2Xuyu9lJZub3\nRzLZyWRmMpNMMpPs83699rXJ5Mw5z3MmOc88Z855Ho8kSSAIgiAIIhWv0wIQBEEQhFshI0kQBEEQ\nOpCRJAiCIAgdyEgSBEEQhA5kJAmCIAhChxKnBcg3gUAwZ8t5a2snoK9vOFfVOwbpVXgUq26kV2FR\nSHr5/TUerePkSdpISYnPaRFyAulVeBSrbqRXYVEMepGRJAiCIAgdyEgSBEEQhA5kJAmCIAhCh3G3\ncIcgiPHN4OAAAoEup8VIoq+vGr29IafFsB236DVpUi2mTJkCj0dzbY4h5EkSBDFueOutN7F9+zan\nxUjBDYYkF7hFr337mrF8+VIMDg5YPpc8SYIgxgXRaBSHDnXgm9/8ttOiEHlm5kxg1qzZWL58KS66\n6BJLHiV5kgRBjAsCgS4cffQMp8UgHMLj8eDjHz8ePT09ls4jI0kQxLhgeHgIVVVVTotBOEhdXR0G\nBvosnUNGkiAIghgXZLBuh4wkQRAEQehBRpIgCMIkvJCz0M+uYM2a53DddVfbUld3dwA///mP0pZ7\n7rlnsXLlU7a0mQvISBIEQZiAFyRcs6W76A2lXdTX+/HPf/47bbmdO99FOBzOg0SZQVtACIIgTLCi\nOYQP+yNY0RzCJcfX2Fr3229vx4MP/gPTpx+BvXubEYlEcM01v8VnPzsLoVAIixffgaYmDh6PB6ed\ndgZ++tPLUVJSgi9+8XR8/vNfwJ49jbjpplvwy1/+BPPmXYyGhk0YGhrCL395FV599RXs3bsH9fV+\n3HHH3aisrMTzz6/CqlUrEY1GMDg4iEsu+QG+853v6sq3Zs1zeOWVFyGKIrq7A6iv9+P3v1+E+no/\nuroO4c9//hM6O9shSRLOOecbuPjiS9HR0Y5LL70AL7+8CQ8//C90dnagp6cbnZ0dmDy5FjfffDve\nf3833njjdWzb9ibKy8sxa9Zs/OlPN2N0lAcg4Rvf+Dbmzv0fW/vaKuRJEgRBpIEXJLzaPgIAeLV9\nJCfe5AcfvI8LL7wEjzyyDF//+rfw738/AAD461/vwsSJk/DYY0/ioYcex549TVi+fCkAIBKJ4Mwz\nz8Ly5StxwgmfAM/zmDKlHo899iS+853v4o47bsFVV12LpUufQigUwqZNGzE8PIznnnsWf/7zPXjk\nkWVYtOh2/OMf96aVb9eu93DNNb/F0qVPgbET8de/3gUAuPnmP+Czn/0cHnvsSdx//7/x4otr8cor\nL6ac/9577+CPf/wTli17GjU1NVi1aiW+8IUv4vOfPxvz5l2M88+fh2XLHsMZZ5yNf/97Ke666x68\n997bEEXRxl62DhlJgrBChHdaAsIBVjSH0BKKAgBaQlGsaLY/ksy0aYfhuOMYAICxExLRYbZubcD5\n58+Dx+NBWVkZvvWt8/Hmmw2J8z7zmZOT6pkz50sAgCOO+BhmzpwJv38qvF4vpk+fjmBwABMmTMCd\nd96NhoY38OCD9+Oxx/6NkZH0OR9POeU0HHXU0QCA8877Dt58cytGRkawa9d7mDt3HgCguroa5577\nDWzd2pBy/sknfw5VVdUAgOOPP0Ez+s3ZZ38Ry5Y9huuv/w02bnwVV1/9G3i9zpopMpJ2QAPn+CDC\nw/eX24FIxGlJiDyi9CJlcuFNlpeXJ72XJCn+X1QdFxGNRhPvKysnJH1eWlqWeO3zpT5R6+o6hB/+\n8Hs4dKgTn/70/8P8+b8wJZ/PN5YbUhRF+HxeSJKYkHPsMylJPhk9/ZSceeZZeOKJlfjSl76MpiYO\nl156IdraDpqSL1eQkcwWGjjHDZ6X1sKzrxmel9c6LQqRR5RepEyuvEktZs8+HStXPgVJksDzPFav\nfgannHJqxvV99NGHmDx5Mi677Mc49dTT0dCwCQAgCILheTt2bEsEhn/22adx5plnYcKEKnzykydh\n5coVAIBQKIR1616wJJ/P54MgxPp34cIbsH79y/jyl7+Ga6/9HaqqqnDoUGcmatoGLdzJkqSB87Lv\nOS0OkSsiPLzbtwIAvNu2QPjKOUBpqcNCEbmGFyS8eFB7KvLFg8OYN7MaZb4Mdqhb4Oqrf427774L\nl156ASKRKE477XRcemn6rRV6zJ59Gl54YRUuuuh8VFZW4MQTP4nJk2vR1nbA8Dy/fyr++Mcb0dPT\njRkzjsF1190AALjxxluwePEdWLPmOUQiEXz1q+fg3HO/ic7ODlPynHbaGbj77jsBAD/4wU9wxx1/\nxKpVK+HzeXH22XNw8smfy1hXO/BoubzFTCAQtE/hCA/fbTfB094GafoRqL37bnT3u3cpc6b4/TUI\nBIJOi2E7VvTyvLAKvtVPJ94L3/oupHPPy5VoWUPXLJXm5iYAwMyZx5k+RxAltA/re1jTJ/jg8+bW\nSLqBNWuew2uvrcedd/7VaVGywug74PfXaF5Imm7NAs9La+Fpb4u9bm9DeNWzDktE5ASFFynj3baF\nptjHAT6vB0dWl+j+jQcDOd4hI5kpGgMnv2kTDZxFiPJmKHGsvY2eTRLjhnPP/WbBe5GZQkYyQ7QG\nTvHAARo4i41IBN4tmzQ/8ja8TjdFBFHk0MKdTEgzcNKijiLC64VwxbWGnxMEUbyQkcwEnYGzrq4K\nvb1DNHAWEz4fMO1wp6UgCMIhyEhmgs7A6fPXAKXFt6KQIAhivEIuD0EQBEHoQEaSIAiCIHQgI0kQ\nBGEWF8VpViYrfvjhf2Hx4jsAAL/+9ZXYt2+vk6IVFWQkCYIgzOCyOM16yYr//Od7ccwxxzogUXFC\nC3cIgiBMoIzTnIuQhKtWrcR//vMEvF4f6urqsGDBdXj88UdwzDEzcfHF3wcA3HrrQhxzzEwcccTH\nkpIVK/nud7+JW265Ayec8Am88cbrWLLkYUSjEVRUVODyy6/GSSd9Gg8//C+8//4u9PR0Y+bM43DZ\nZT92XbJjt0BGkiAIIh05DnC/Y8c2LFv2GP75z0dQW1uLNWuew/XX/xonnPAJzfJf+MIX8cYbG3HM\nMTNx/vnz8PDD/0opc+DAfjzwwH3429/+hUmTJmPv3mYsWPBLPPFELHxmZ2cHHnvsSZSUlOD222/G\nGWecje9//wfo6enGvff+Bd/+9vmO53J0A2QkCYIg0qCO02y3N/nmmw340pe+gtraWgCxMHD33PNn\nXSNphm3b3kRPTzeuuuqXiWMejxcHD8ayfXzyk59CSUnMBJx99hdxyy034cMP38esWbNdkezYLVAv\nEARBGJGHAPeimJqcSJIQz6U49lnEQpuiKOBzn5uNRx9dlvj7178ewbHHzgQAVFZWJsq6MdmxWyAj\nSRAEYUA+Atyfeupp2LDhZfT19QEAXnhhNSZNmoQTT/wEPvroQwBAf38/du58N3GOMlmxFp/97Cl4\n662taG1tAQBs2fIGLrvsIvB86gpdNyY7dgs03UoQBKFHnuI0n3LKaZg372JcddXPIYoSJk+ejDvu\nuBtVVVW4+eY/4KKL5uLww6fj5JM/mzhHmaxYi2OPnYnrrrsBN910PSRJgs/nwx13LE7yIGXcmOzY\nLVDSZRuhRLeFRbHqBRSvbvlOugxBALq79D+vnxoLU0kUBJkkXSZPkiAIQg8KcD/uoWeSBEEQBKED\nGUmCIAiC0IGMJEEQBEHoQEaSIAiCIHQgI0kQBEEQOpCRJAiCIAgdyEgSBOGqPIluRhBEp0VI0N7e\nhhtu+I3TYmRMoeTDzOs+ScbYqQDu4DhuDmPs/wH4J4AogEYAP+E4TmSMzQfws/jxWziOe54xVg9g\nGYBKAO0Afshx3LBW2XzqQxBFQTxPonDt9bZmtig2BEHEW68049Qvz4TX57x/0dnZgf37W50WI2N2\n7nwXxxwzM+X4n/98rwPS6JM3I8kYuw7A9wEMxQ/dBOBmjuPWMMb+D8DXGWPbAFwJYBaACgBvMMZe\nBnAjgGUcxz3KGPsdgJ8xxpZrleU4bjRfOhFEMZDrPInFwr4PAxjoGcHeDwP4+EnTbK9fL/fjbbct\nwsjICP74xz9h795mXHnlz3Hvvf/EHXfcgkAggGuu+RV+85vrcfnl83H00TPQ0dGBv//9AXR0tOH+\n+/+GcHgEHo8XP/rRT3HmmWdhzZrn8NprGzA6OorOznZMm3YY5s79Hzz99AocOLAfF1zwPVx00SUA\ngOeffxYrV/4HkiRi4sTJuOaa63D00TNSZNfKhXnUUUcn8l8Wcj7MfHqSzQDmAng8/v4dAHWMMQ+A\nGgARALMBbI4bulHG2B4AnwbweQC3xc9bG3/drFN2W570IYjCJ8d5EosFQRDR0doPAOho7cexJ/pt\n9SaNcj8uWHAdfvSj72Ht2uexbNljuPLKa3HssTPx29/+HnfffScWL/47Ojra0dV1CDfddAs+85mT\nMTg4iNtuW4TFi/+Oww+fju7uAH760x8kwrHt3PkOlix5An7/VFx66YV45ZWXcM8996O5eQ9+9rMf\n4oILLsZ7772DtWtfwD/+8RAqKirw1ltbccMNv8HSpU8lya6XC/Pxx1fo6ltI+TDzZiQ5jnuaMTZD\ncagJwH0Afg9gAMBrAL4bfy0TBDAJwETFca1jyuOG1NZOQElJ7mIt+v01OavbSUivwsOMbiNPPYWw\nIk9idcN6VH7X3RnpM71mfX3V6O0NZXTuvg8DGBqITVINDYza7k0a5X487rjjsWjRbfjpT3+Ar33t\nXHz1q/+tWYfP58MnP/kpAMD77+9ET08P/vd/f51URo5desIJn8C0aYcBAKZPn47Zs0+D1+vFEUd8\nDDw/inA4jC1b3sDBgwfw85//KHH+4OAgBgcHMHHi2FCrlwuzo6M9Z30CZJ4Ps66u2tJ3yMnYrfcA\nOIvjuPcZY5cD+AuAFxHzKmVqAPQDGIy/HtE4pi5rSF/fsC3Ca0FBpQuLYtULMKlbhIdv40YoozqP\nvLYRoTP+y7XeZDbXLFMDqfQiZez2JuXcjzfffHvi2KFDnaiv9wMA9u9vxaRJk9DUxCESiaBU4/qU\nlZUljIYgiDj66Bl48MElic+7uwOYPLkWL720FmVlZUnnyucpEQQRX/vaufjlL6+MyyiiuzuAmpqJ\nKtmllBjvkgREo1F4PB5kmw9Tq09ef/1VzXyY27a9iR07tuGRRx7EP//5bxxxxMdS6u3tDWl+h/QM\np5NPn3sRM3RAbDFOLYC3AJzFGKtgjE0CcCKA3QA2Azg3XvYcAJsMyhIEYYJ85EksBpRepIzsTdqF\nUe7Hjo523HPPX3D33ffhqKNm4P77YwtbfL4SRKPa+SQ/+clP4eDBA3j33bcBAE1NHC688Dvo7jYv\n8+zZp+GVV15Ed3c3AODZZ5/GVVf9IqWcXi7Mj33sSEyePLng82E66Un+BMATjLEoAB7AfI7jOhlj\n9yJmBL0AbuA4LswYuwXAkvhq1m4AF3McN6RV1hlVCKLAyFOexEJHFES07e3T/Kxtb59t3qRe7sfS\n0lIsXHgDLrro+zj22I/jmmt+i8suuxCzZp2Kk076NHw+H+bPvxSLFt2eVF9tbS1uvfVO3HffPeB5\nHpIk4g9/uBmHHWY+o8mpp56O733vMixY8Et4vV5MmFCFW2+9K+4djqGXC9Pr9eL88y8o+HyYlE/S\nRop1+o70KjzS6lageRLznU9SFCWMhPT3kFZWl8Hr1UxDSLgQyidJEIQ5KE+iKbxeD6omlqcvSBQt\nzu+IJQiCIAiXQkaSKBwodBpBEFmQydNFMpJEYRAPnQYLS8gJQsmkSbUIBOxbkUoUHvv3t8Dvn2rp\nHHomSRQEFDqNyJb6+nq8+eYWbN3agPp6Pzy03mZcIEnA0FAI+/e3Ytq0w5MCIZiBjCThfih0GmET\nX//6N9Hd3Y2BAe1tHU5RV5d5NCA34wa9PB5g6tRpSRF6rEBGknA9yk3v8mZ38iaJTKmvr0d9fb3T\nYiRRrNuRikEveiZJuBuFFynj3baFnk0SBJEXyEgSroZCpxEE4SRkJAn3kiZ0GnmTBEHkGnomSbgX\nrxfCFdcafk4QBJFLyEgS7oVCpxEE4TB0K04QBEEQOpCRJAiCIAgdyEgSBEEQhA5kJAmCIAhCBzKS\nBEEQBKEDGUmCIAiC0IGMJEEQBEHoQEayiOGFDDKM5plcypg3/SkZNEEULWQkXUq2AzwvSLhmS7er\nDWUuZcyb/pQMmiCKGjKSLsSOAX5Fcwgf9kewotm9OepyKWO+9FcmgyYIovggI+lCsh3geUHCq+0j\nAIBX20dc6U3mUsa86a9KBk3eJOFWhHDYaREKFjKSLsOOAX5FcwgtoSgAoCUUdaU3mUsZ86W/VjJo\ngnAbQjiMbcu3QAyPOi1KQUJG0mVkO8ArjayM27zJXMqYN/0pGTRRILS8sBX95VOw74Wt6QsTKZCR\ndBF2DPBKIyvjNm8ylzLmS39KBk0UAkI4jI7B2DDfMeghbzIDyEi6iGwHeF6Q8OLBYc3PXjw47Apv\nMpcy5k1/SgZNFAgtL2xFqHwyACBUPpm8yQygfJIuId0AP29mNcp8HsM6fB7gttlTDD93mlzKmDf9\nKRk0UQAkvMjysWMdgx4cEx6Ft6Jc/0QiCTKSLsGOAd7n9eDIandf0lzKmDf9KRk0UQDEvMi6pGOy\nNznz/C84JFXh4e4RdRxRCAaOIIjCQAyPom2wJMmLlGkb9JE3aQEalQmCIIoNnw+fm3O0zoe1sdkQ\nwhRkJAmiWInwQGmZ01IQDuAtLUHVMUc6LUZRQCsMCKIIkXiKKUsQdkBGkiCKkPCqVRRTliBsgIwk\nQRQbER78G7F9nBQFiCCyg4wkQRQZnpfWQjxwIPaaogARRFaQkXQJboiGQxQBFFOWIGyFjKQLKIQE\nyURhQDFlCcJeyEi6gEJIkEwUABRTliBsh/ZJOow6f6SZGK0EoYkipmxdXRV6e4dSPicIwhpkJB1G\nK3/kJcfXOCwVUZAoYsr6/DVAadBhgQii8KFbSwcphATJBEEQ4xkykg5SCAmSCYIgxjNkJB2iEBIk\nEwRBjHfomaRDFEKCZCJzeEGiBVgEUQSQkXQIyh9ZvMj7XhefXk+GkiAKHJpuJQiboX2vBFE8kJEk\nCBtR73vN+bPlCJ/b+glinENGkgCQeexYpxcYqds3K4+dcivr0tr3mjMilDNSl3zePNCNSlGTVyPJ\nGDuVMfZa/PVUxtgqxtjrjLHNjLGZ8ePzGWPbGWNbGWPfiB+rZ4y9xBjbxBh7kjE2Qa8sYZ1MY8c6\nHXNW3b5ZecbKibbKkO99r56X1lLOSC3yefNANypFT96MJGPsOgAPAaiIH7oTwP9xHHc2gN8DOIEx\ndhiAKwGcCeBrAG5njJUDuBHAMo7jzgLwDoCfGZQlLJLpMzSnn72p2zcrj1xuye4+W2XI675XRbYP\nyvKRTD5vHuhGpfjJ5/LKZgBzATwef38mgJ2MsVcAtAC4CsB/AdjMcdwogFHG2B4AnwbweQC3xc9b\nG3/drFN2m5EQtbUTUFLis1OvJPz+wgopNxoV8XpXNwDg9a5RXH7aYSjzpd47qfUye16uULf/k1Om\nmpJHed6avYO47KSjM5ZbWdfGQ6O6nukrHWHb+2fkqacQjmf78LS3obphPSq/+z9JZQrtu2gWI70k\nnsfgO29BBFD69puYePE8eEpLcyKH3W2Nx+tVCOTNSHIc9zRjbIbi0AwAfRzHfZkxdiOA3wJoBDCg\nKBMEMAnARMVxrWPK44b09Wlv4LcDv78GgUBhxctc2hhEc3/smUpzP4/7tnamxI7V0svMeblE3f6v\n1h0wJY+dcivr2jvAY+6MKtwyq06zbG93CD6vTdtBIjx8GzdCWdvIaxsROuO/gPggXYjfRTOk08vz\nwir44gmnxQMH0LNsBaRzz8uJLHa2NV6vl5vQM+ZOLtzpAbA6/vo5ALMADAJQSloDoF91XOuY8jhh\nkkyfoTkdc1ar/V29yYsntOSxU26tut7uHsW0Sh+OrC5J+bPNQIJyRuqSz4TTlNx63OCkkXwDwLnx\n12cDeB/AWwDOYoxVMMYmATgRwG4AmxVlzwGwyaAsYZJMn6E5HXNWq/2oys5pyWOn3I71AeWM1CWf\nNw90ozJ+cDLky7UAHmKM/QKxadOLOY7rY4zdi5gR9AK4geO4MGPsFgBLGGPzAXTHyw5plXVGlcIj\nXexYvbyWmZ5nF0btG8ljp9yO9oEiZ6Te5+OSNDcPwlfOSUxFF1RbhON4JGl8BdIOBII5U7iQ5t8F\nUUL7sKD7+fQJvsQUoVIvK+flAnX7gighEB57P7XCB6+ifVkeLbmViYmtyO10H5ihkL6LVtDVSxCA\n7i79E+unxvJt2kEO2hp318uF+P01mj9aCh46Tsk0dqzTMWe12p8xMf1du9Z5/kllqOJHbZGBcBhF\nwumiaotwnHE6N0MQBEEQ6SEjSbgGp0PcmaVQ5CQIInvISBKuwOkQd2YpFDkJgrAHMpKEK3A6xJ1Z\nCkVOgiDsgYwk4Th5Ty+VIYUiJ0EQ9kFGknCcvKaXyoJCkZMgCPsgI0k4itMh7sxSKHISBGEvZCQJ\nR3E6xJ1ZCkVOgiDshYwk4Rjpwru5xUsrFDkJgrAfChtCOIbPA9w2e4rh526gUOQkCMJ+yEgSjlEo\n4d0KRU6CIOyHplsJgiAIQgcykgRBEAShAxlJgiAIgtCBjCRBEARB6EBGkiAIgiB0ICNJEARBEDqQ\nkSQIgiAIHchIupxcR3Oxo36jOgoxGo1SZvm1WT2y1dctfZkvfbMmwjvbvhq3yUNkDRlJF5PrBL92\n1G9URyEmKFbKLL8O8aIpPazqqy7nlr4025bj1zfCw/eX24FIxJn21bhNHsIWyEi6mFwn+LWjfqM6\nCjFBsVJm+fWiHb2m9LCir5aBcUtfmm3L6evreWktPPua4Xl5rSPtq3GbPIQ9kJF0KblO8GtH/UZ1\nuD1BsZYXp5R5fdswNsRf7+yNTaEZ6RHiRUv6qg2MW/rSbFuOX98ID+/2rQAA77YtzntvbpOHsA0y\nki4l1wl+7ajfqA63JCg2nroUk94v3xNMyLx/SEBr/HU0XoWeHrwg4Scbu0zrq2Vg3NKXZtty+vp6\nXloLT3tb7HV7m+Pem9vkIeyDjKQLyXWCXzvqN6rDLQmK9Z6ZyV7ckt19Se9Xt2qnw1KipcfyPUEE\nRsW05ZTtKw3MsqagK/rSbFuOX1+F1ybjqPfmNnkIWyEj6UJyneDXjvqN6nBLgmKtZ2bKAX7N3sGk\nadJ+XtSsR4laD16QNI2rkdepNjCrW4dc0Zdm23L6+iq9tsQxB703t8lD2AsZSZeR6wS/dtRvVMe6\nA0NY54IExXrPzJQDfHM/j4U7elMG/HQo9Vi+J6hrXLX01TIwAxHtPslnX5r9XjidgFqKRODdsknz\nM2/D6/n33twmD2E7lCTPZeQ6wa8d9avr4AUJf3q3D/97ci28ADwAvF7tivKVoFjrmdm8mdUpXtyu\n3uR9bdMrvfjJiRNR7vWgNz6FOqXcm9DHX+GDz+uBzxPTe2NHOPn8CT4smlWH0nh5pb5GBmZqhRe3\nzp6SOA8ABFHKW1+a/V44noDa64VwxbWGn+cVt8lD2A4ZSZeR6wS/dtSvrmNpYxB7g1E0dIZxyfE1\n2YqYNXrPzARRSvHioirHp31ExP3vD+KxL01DWZoRf3ljMLG4J3H+sKDbD+kMzPQJMQPsBGa/F04n\noPb4fMC0wx1rPwW3yUPYDt3mEFnh+FYADfSemT3TMmTq/MCoiCf2BA3LZDLtKBsY+W9apS/pvVMG\nkiAIfciTJLJCa1oz394kL0gJr8/IeFWWeLD4jHqUej0QRAl8eTkGBmJlp1b4EJWAhTt60T4s4LWO\nMC78eI2uN5nttKO88nbx6fVpPVaCIJyDjCSRMXrTmvNmVudt4FcbGytTmn5/FQLlY4tuljYG0T4s\nAABa0xj8bKcdlStv3TBFTRCENjTdSmSM01sBZBmU2zzUU5rqP70pTaf2I7plipogCG3ISBIZ4fRW\nAFkGu4yNU/sRnYxGRBBEemi6lcgIx7cCwL7noekMvp3Tx45OUUd4oLQst20QRJFBRpLICKe3Athp\nbJQGf3XLEFa2DGHujCqcN6Mq8bldGHmsOX02GU/jJFx7PVBamrt2CKLIoOlWoiCxc3pUNvjTKn14\nu3sUAPB292hii4ZdWzOcnKKmNE4EkRnkSRIFR66mR3O9ncWxKWpVGifhK+eQN0kQJiEjSRQcuTA2\n+XhW6NQUtVYaJ+nc8/IuB0EUIjTdShQcmW7zMMIN21lyAqVxIoisICNJjHvcsJ0lV1AaJ4LIDppu\nJcY9btjOkhPSpHGiZ5MEkR4yksS4x+ntLDmD0jgRRNYU4chAEAQASuNEEDZAt5IEQRAEoQMZSYIg\nCILQgYwkQRAEQehARrJAsbotIV15u7Y5qOvJtF475JHrMKpL+VmmspstF+LFpP96dcl/ZuqW67Kj\nn5XtZlKPXXVlSyFv2SHcBxnJAkRONGxlEDcqb7U+s+1YqVc9wGYjjzxAX7OlGyFe1K1rNDr2Waay\nmy0X4kVcsL4Th4aiuGB9p6ah5AUJVzcEcNXmAK5uCBjKrqyzLyxk1F9K2eW2FzRkVs+Chm5c3RBI\n1KV8n0/s+i67hgjvtATjnrwaScbYqYyx11THLmaMbVG8n88Y284Y28oY+0b8WD1j7CXG2CbG2JOM\nsQl6ZccD6kTD2Za3Wp/ZdszWqx7Y5POWNQUtyyDXtXxPEB/2R7BoR6+uDEt29yU+MyO71sBrVseF\nO3oRFoD5m7oQFoBFO3o16+IGomgcjIIbiGKhgezKOq/YHMjo+illl9v+aCCzej4aiIAbiCbqUr7P\nJ3Z9l11BPHMLRUdylrwZScbYdQAeAlChOHYygB8D8MTfHwbgSgBnAvgagNsZY+UAbgSwjOO4swC8\nA+BnBmV5M3uhAAAgAElEQVSLGquJhtOVtytxsbqeEC+arlc5sCnreXJvyHBq0qiu1a2xCDo7e/kk\nGZSe4tp9gwCADe0j2NA2bCi7loditu9CvIhdcTmG4pHvdvbySbrxgoT1bclRf3apZNers2NETCuD\nGqXsG9pH8MrBocRnGyzWs0ER8/aVg0NJemxoG87b1Ktap0L3JilzizvI5z7JZgBzATwOAIyxKQBu\nA3A1gAfjZWYD2Mxx3CiAUcbYHgCfBvD5eFkAWBt/3axTdpuRELW1E1BS4rNTryT8/hzmBATwwLs9\nSZkq1hyKYP5n9KPFpCtvtr50eqnruW3XgKl6R6MiXu/qBgC83jWKisqyxHm8CNy+awAP/PdRhm1r\n1dUfN0DR+DjZEopidQePN9qG8Mg5R+LRXX1o7o8ZmVZFzFY92SUJ+LA/kqSH2b67fu3+hBwyUSlZ\ntwfe7cH+ISGljF7dWnWqyxldM6XsraqYta0mvlfKepTnHxhOvqlpHRKwun0Ub7QP45FzjkSZL/v7\ncj291Dqt7hjFL072Z91evlDqJfE8Bt95CyKA0rffxMSL58FToNGRcj0m5pq8GUmO455mjM0AAMaY\nD8DDAK4BoEy9MBHAgOJ9EMAk1XGtY8rjhvT1acfotAO/vwaBgPUpQrPwgoTnm/qTjj3X2I9zp5Vq\nZqpIV95sfen00qpnR2dyRg09OZc2BhPGqrmfx7Kh5GnI7Z0j2Nc2gOqy9IOrsi4tlr3fi4GIhL9t\n6cDGjrBuObXsqxv7AUlK0gOAqb4L8SJ2HEquT0bWrcznwerGPkPdlHUb1SmXO+KwibrXTOt6qVlt\n8L1S1vNcmnoAYNkHfRiISLhva2fWqcf0votaOi3/oB/nHV5uWxaXXKLWy/PCKvgOHAAAiAcOoGfZ\nioLM3JLrMdFO9Iy5Uwt3PgfgOAD3A3gCwCcYY38FMAhAKWkNgH7Vca1jyuNFi9VMFenK25X5Qqse\nLS9HXa9WeqqBSPKJUUn7+Z0arbrUyHWvah1OkVfdppLWUBStcS9P1sNs3y3c0ZtSn7KdRTt6saI5\nlOJFqlHWbVSnmeunJbuaVpP1qL1QLeR+z2Y6Px1aOvXzIp7YUxgDdBKUucVVOGIkOY57i+O4T3Ic\nNwfAhQA+4DjuagBvATiLMVbBGJsE4EQAuwFsBnBu/PRzAGwyKFuUWM1Uka58iBdtyXxh1E66es0M\n1kDq8zstzNYFjE3FZsqG9hGsOzCk+ZlSxxAvJp6J6vFeL481reYG8hcPDqMvLKStMyaDto5Wrtc6\ng+8BL0hYZ7IemVylHjPSaWXLUME9m6TMLe7CVbFbOY7rZIzdi5gR9AK4geO4MGPsFgBLGGPzAXQD\nuJjjuCGtso4Jn2OsZqpIV77Ma0/mC3U7oiihKzzmFfkrfEn5HeV6rQzWUQl4piWE7x8/UfNzo7rq\nyjy44lOTcf/7A+gKjxmOmhIPggp37CtHVGLO9ApAQkLeqRU+PL9/GCtbkg1iayiKuTOqcN6MKs02\nZR3LvMCXplfi5bYxD/fkulL8v/rY+rJja0rg83hQW+FFb1hE76gIUZKSjHhtmRdTKnw4fEIJvF4P\nKrzArafUQRQl9IymGkI2uRRlPi+8Hu0LqLxeoiihYziKXkU9dRVe+Dxj+ut9D3we4NZZdYlrLYgS\nusNCQvbJZV7Ulvvwd1W/253IWqlTRJSwcHsv2obHvn8TS70QRKlwUrlQ5hbX4ZGkwrrLypZAIJgz\nhQtp/t0KudBLECW0DydPMUZECTe81ZM0qMpMn+DDw1+Yqjm4atWl5NW2YSxpMvZgZlSX4P6z/En1\n84KEH2/s0qzbSB7l+b/YFEjycLXayQXZXjNekLKWcWljEI80psrww+NrMn42aaRXLtrLFwm9BAHo\n7tIvWD81Fri+QCikMdHvr9H8wrvKkyTGD1rpqQRRwp2n1eufo+fVGKS64gUpyZPTQ54KVA6m2eaZ\nNHpu6eZBW97usvj0+oShtGo000332+1N5ru9nEGZW1wHGUnCNeQir6Pa0Bl5q+rBNBt5CnnQVu5b\nveT4Gk2jmY58J7Iu2sTZhOOQkSSKGrWhU3qr1ZMqceXLB9E+LGD6BB8WzapLGkyzmXIs1EFbHSBh\n3szqFKNphnwnsi7axNmE41DsVmJcIQ+mR1aX4LX9Q4nnje3DAho6w4lFO9nGAFW2o/WnXMzkJpRT\nxC2hKJY1BW2JyEQQhQoZyXGG1TBvZijEgVMZlk5GaQSKKgaoSbT2mq5uHUoymuOpPwgCICM5rpCz\nRthpKAs168KK5lBKdB7ZCNgVz1aN2/tIa6GROrjDeOoPggDISI4r5KwRZiLYmKUQPa50i2qWNwVt\n957cfjNhdt/qeOkPgpChJ93jBGXWCDmCjZl4qEZoLfJw64pNJfKimrq6KvT2jgUKiIixAfsW1U2E\nHbplsvgln6gXGg1HRVy5uVsz/J0dq3Pd3h8EIUNGcpygjPcpxwy963T9PYlmUC/yKJQBT15U459U\nhip+FMCYZ3OKvzwRp1UmW90K4WZCvTr0MW4QUQm6kYWyEb8Q+oMgZGi6dRyg9CJlzMRDNUJrkUch\nr36UPZtnWtLHZM2k7kJa/MILUiJLytvdo5hW6bN1dW6h9QcxviFPchyglTUiW2+yUKPJaKE0+BNL\nvVh8eh1KNYxAJs6O3s2Em72nXM4QFGJ/EOMbS54kY+wzjLHP5koYwn6MMlG8l6E3aTUjidtZplio\n0xbfL2nX3ka70pHli1zPEBRafxCEaU8ynonjawCmADiWMXY8gMM4jns9V8IR2VMWzxph9LlVCjWa\njBYhXsSTe5MHaLs8GzOh6QC4yoPK5QxBIYfqI8YvVqZbzwdwEoDt8fcHATwK4AybZSJspKzEi1Om\nVthaZzGFAFu4oxdqZ9ouo5DuZkIQJfzmzR5LMVFzSa6NWDHdXBHjBysjXQ/HcQJjDADAcdwwY8ze\n0Zcg8kiQj6YsaJKxxSikuZlY2hh01TaIXBuxYrq5IsYPVr6xuxljXwIgAQBjrBJAdU6kIog88H8f\n9CctaFJvd8ilZ+PGbRBkxAgiFStPpH4H4DoAM+LPJ18D8HIuhCKIXMMLEl7cl5wMVr3dIZdByGkb\nBEEUBqaNJMdx/RzH/TeACwAMArgHwBW5EowgcolR7NZcU2x7TAmimLGyunU5gAUcx70M8iCJAsbp\nVZZ2rSBV5rvMJvelGyk2fYjCxcoDiJM4juuU3zDGDgPwV47jLrRfLMIJ1ANTsQ5UerFbgZjOsspa\n+qsNk4xRP8nnyOWtGGj5HHX9IV7E796KrYzlBSnxWqtcdZk3Re4QL6Ku0qepl57uRnqlK2sUK1hd\njxwi8E/xRURGMYbTXSOCyBYrzyST5qbiBvN4e8UhnEKdlaFYsjTwgpSig7xAZcaksqRgAdMqfbh7\nVz8ESVt/5bEQL2JBQzd+9cYhLGjQ7yf5nBAv4pot3RBECbfNnoJH50zFo3Om4r4z61HuBf7x+Xrc\nNntK0mIhXpCwoKEbVzcEUuT48cYufNgfwRN7gvhJ/LV6qjjEi5j3Sid6R4SE3Lwg4Yo3ujBv/SH0\nhYUUvbTep9NLLZ8avRRtsjzK/uEFKREi8IevHTJM7ZbuGhGEHVgxkhsZY4vkN4wxH4DUyMdEQaJO\neVWIKbDU6BkZmdFo8uCr1FlLf/nYE3uC+PHGLnw0EEFzUMRHA/r9JJ+zaEcvPuyP4Ol9Q0mG+aGP\nBjEqAg99OJiyWGhFcwgfDUTADSQ/K13WFET3aEz2Z1uHEIi/Vj/XvHF7D0ZF4IrNgSS99gQFSIgd\nV+ut9d5Irxu394AbiGJ5U1CzHC9IminaZGO2fE8wqX+WNQUTz2t7eckwtZvRNSrk7y3hLqwYyT8A\nOJ0x9hZj7C8A1gN4IzdiEflEvR0hxIuqpMP2JWnOJ3pGBojp/KN1B5K8J1nn9W3D2KBKuqz8/JmW\nUMJIyWzQWHijPEcODag0ZFrpy5TnblAs7tnQNpyQY3Xr2BTxAD/WpnLh0eBoFDt7IwCAzrCY0Ovl\ng2PndoyI6BqOprn2qTcXSr12xdtY1Tqkqf9VmwPY2ZOqo2zMVrcOJ/XP6tahlOe1WsH4tRJj5ypZ\nNjG+sbK6dYjjuK8CuAHAIQAPAfhZrgQj8od6O8KiHb1J75fs7nNSvIzQMzIyK5pD2BUIJ3lPss77\nhwS0qrZnKD8fjKS216qxMlZ5jrwfU2nItNKXKc9tVRiL1iEBK5pDWNYUxEBEf/CXjcOC9W1Ql9o/\nJODgcLKxmf96l+G1V+ok959SL7m2gYiU8CaV5RoHo5ATj8k6Ko1Zf9z4yf2gpZu6b9QyaF0jMyuV\nyYgSZkhrJBljL8T/L2CMfQHAmxzH3clx3FKO46JpTidcjtZ2BHVA9DV7BwtuQNEzMoCx56yF0rM0\nQulNavWrzKvtI+gdEXTTl6kNvMwrB4eSvEgtWkJRPM4NYsehcFp5AUDltKVce6WXJj871NNrVetQ\n4tliiBexvi11gdJ7PTyWNg6meIvpUHqTWn27oX0EG1TtGXmT9OySMIsZT/Ky+P8jEPMi9zDGmhhj\nTzLGfps70Yh8oLUdQZ1Wq7mfL6hnPHpGRvYmjTxnLZSepRFKb1KrX2VaQlFc0RDQTV+mNvAyB4ZF\nQy9S5qm9Q8h0glwtk9JLk58d6uk1EJFw07YefNgfwcIdvdivSl4NAEJcvkzkkr1Jrb5tDUV1k2Vr\nQc8uCbOkNZIcx3UzxjwAzuE47qscx00F8F8AloMW7hQ0RvsF1RRKCizZCGoZmdYhAU/sCab1nLNh\n3cFhhHgxbb92jmibsfd6eazZn/nAXV/ugcZscFasOzCUuOlI11c7+2Kt68XEBQA+w6/Re708ekci\npr+zgPb3lp5dElbwSJK5LwhjbB2Ab3EcN5pbkXJLIBDM2S/C769BIKC9ys+NCKKE9uGxu29RlNAV\nHnvvr/DB5/Uk9hNOn+DLaai2bOEFCVc3BDDIi+jQMULVJalTjEAsbuvXj5qAQDjZG5lS7kXnSOzY\n7e/0YTjVOUpQV+bB7bOn4KjqEhwKiyn9KTOpzIM9A1FMKffCq+pPUZRQV+FFv8qSCKKE3lERNaUe\n/OvDQXSMiJha4cUVJ01CqdeD2jIvyku8GI4I+P22PvRqbJso9wJfPaISpT4Pqko8GIq7jRNLPThl\najmCESAiSPjb+wPoCouYPsGHRbPqsHb/MFa2jHl/3zqqEps6w+i1aO08QMpzUmXfXfmpSegLCxjg\nRXg8MZ08nrH+mVLuxRw2BY1tgyl9Jvez/J1Vov7eLm0M4pHGsd/pD4+vcTzAfKGNHWYpJL38/hrN\nwc1KMIGDAFYyxq7gOG6vPWIRTqIV0ProiaUp5fyTylDFu//eaEVzCNxAFN8+egJOm5aaoCYiSrh3\nVz9C6jlFAFu7wph/4kTMUOnPCxIW7xrAbafUobrUi2Gdlb6HV3px8ylTcHR8G8eR1bFJGq3+BIDj\nJ5dbVQ9AbICXbwC6wiL2DkYxb2Y1rtnSjcWn18M3wYfFZ9RrBkoAUg2GVv1d8dWw7cMCXm8fwdvd\nydf+3d4I7jitHqXe5OADC3f0Jt10qVH2+icmleK6k2styQYA1WU+zSDsev2sRi8koBsCzBPuxIqR\n7AZwFIA3GWMhxPJKbuc47o6cSEYQFlAOfu/28PjZJyalDHqCKOGu0/2J92pDojVGys+uVrUO4c7T\n6pM+W90ylPCwOkZENHSGcWwOPRK9AV6QpKSUW0dWl8A/qQylI2FLA79W/atbh1Keg7aGomjoDCd5\nX0sbg5oGcu6MKnz5Y5W4cnN30vPOxsEIasu8htF0ckEuk0oTxYmVLSC/iz+T9AP4AoD/Az2TJFyC\nmeX/sucs/6kj7qi9GKXReK09nJQhpLbMm+JhmX2+lekzML0BXt5rqGx/NCqaipqTrn69hULKZ31G\nz7a3doXxwPsDmouUFm7v0ZUtF6SL2UvPJgktrAQ4LwcwD8AUxDzIZwE8myvBCMIsuZpC0zK8lxxf\nA16Q8JONXYlINzJmPBJ564FWnFUjjAZ4ea+hsv0lu/sSEWx+cMLEtHIY1T+1wouFs+pwy9t9aB8W\nEs8qZfGNkjUPR0Vcsblb87N3eiPoCwuorfBpfm43uU4qTRQnVqZbVwKoAbAbwCWMsQEA53Mc158T\nyQjCJLmYQjMyvMv3BFMMpEy6LCLKrQdWZNMa4LWeA77aPoJvz6jC2n2xxS1P7g3hu8dWp0xrquVI\nZ0BebRtOtNM+LCRNLRsla+ajIm49pS7l+Pq2EbzcNoJVLSH84IRJJnogeyipNJEJVr4xHwdwAsdx\nEgAwxn4K4B8ALs6FYARhhnRTaJl6k3qG94k9QWzsiG3Ulz2qUtU0rV5z6q0HStnSZa7weT2YVulL\nKr+iOZTyHFDe9ynnyuTF2P7Cu04fe56qJ4euoRMkvNaRHJxgg1lP3ePBKVOTF1HxgoR/fhAz4hs7\nwrj4uImuyD5D2UNyQ6H3q5Wn5gEAiSV5HMc9AODTtktEEBaQPSA5q4byT51VwyxGhndly1BiD6bs\nUakzieit0NR7bmol64Yc/ebqhgDW6cio3suojn1qNXyb1r5TrTB8RjKr61OGAFQGR9fKSJIPKAJP\nbjDzbNztWDGSgwBWM8aOBwDG2FEAKFwF4SjqxTjqv0z2deoZ3gfP9mNiafJPRrlYxmig1Zu+VaaG\nMpN1Q45+ww1EMbu+PEXGuTOqdCP5pJNDC16QdI3xujSLXbT00mpfGRxdeU4+DRdF4MkN8rPxQu5X\nK0byHcTiGW+OP4/cA+BtxtgcxhitnSaKBj3D29AZRpvG9KaZFFN607fK1FBmsm4oY5SuPTiM2jJv\nkhe7tUs7ZqvsTRo9v9XsCw8wu157T+fpUytMTy1rBUeXkYOjq89Z3hTMywBLEXhyAy9IiWfjhdyv\nVraA3MBx3H/Ht4B8GsBFAPoBXA+gOUfyEYQrSPfs0yjFlNG5z7YMpZ36VBoWZYxS+XmjjOwBz52R\nujMrKgH/2Ru0vAViJCphzQH97R16457WlK5RPzzTMoTle4JJ56yKB3PP9QBrdfqZMMeK5lDi2Xgh\n92tGS704jmsF0ArgaXvFIQh3km7158p92ttFjM6VV6cGo8mrU9ULeowylMgeYnWZN7G4R8+bXN8e\n1lxopNRRSy/1Vsm5M6pwXtwQa52jN6V7/jFVWPi5Oty0vSclbGB1qQevqc6R92jmcrM/ReDJDcXU\nr7QemiBMYLjNQZASK15llAOC1rlGq1OVBsEomwgw9rxRXr0qG+RMw9KpZVSvagWAt7tHMf/EiYbb\nXLSmdJ/eNxadSI1erF0ZuT/thiLw5IZi6tf8xoQiiCLE6nM+XpCwoKEb6w5op4ySpz7NZmlRrl6V\nDbI6mlAmC5n0sqmk001P5nUHhjRzTJohF9N1FIEnNxRbv5InSRBZkMk+zRXNIXw0EEmatlQjnyJP\n08qZLuRN+EqiEvBMSwjfPz41sk6m6OXklNHTzWhaWhnrVokyIH1tmRd/2N6bCLKubvNynQDzmUAR\neHKDsl+1ZjQKrV/JSBJEFlgdaJXPatJNWwJImqY9vKoE9+4e0Cz30sERXDCzxrbnPXpeZLrnkXrT\n0rwg6T4rfSswmghIL4hSSiB5JV6PfSMsReDJDcp+LZQMQkbQN4QgssDqQKsXD9ZUW3nyfNIFLE9n\n2LUwK3u6/ixxcT5TojghI0kQeSLbFX/58nxyYYzJayMKFfrWEkSeKJQVf2TQzFHoMUkJc9DqVoLI\nA8W24m+8Q7Fexw90u0jkFTfcfVuRwWxZvXLyca0pzIgoJTb2K0+12kdyMAH5XAAYHI1aqktZh/q8\nEC+CFyTUVabmfZTLKf8DSDmm/K8lV7r3St3M1qPVrhFaOujVo5fyzEg2o3bVr7V0N3NMTycic/Lq\nSTLGTmWMvRZ//f8YY5sYY68xxl5kjE2LH5/PGNvOGNvKGPtG/Fg9Y+ylePknGWMT9MoS7sUNd9/p\nZFAPWGbk1SunPK6OBzut0ofFO/sxrdKXtH/Rah+FeBEXrO9MGLKrGwK4cnMAX36yGb0jgqm6lHWo\nZQjxIr77UifmrT+EvnBy4AO5XIgXE/+vbghgQUPyMeV/ef+nOtOH0XvlsSAfNVWPWjaz1zDEi1jQ\n0I2rGwK69ejFelX3m9l2Y22NXb8FDdq6W/1uuuH3VgzkzUgyxq4D8BAAObncPQCu4DhuDmIJnX/L\nGDsMwJUAzgTwNQC3M8bKAdwIYBnHcWchFmj9ZwZlCZfihkwLRjKoBxWz8uqVSxfw3Oo5Wizc0Yuw\nEIu6I2cHaRqMIiwAVzYETNWlrEMtw8IdvYgAkABcsTmgKeuiHb2J/9xAFB8NJB/7sD+ChfH/chYT\npVzp3iuPXbOh3VQ98mtlu0Yoy380EAE3ENWtRy/Wq1x2WVMQP9nYZardZU1BcANR/Hvn2PX7aEBb\n93TH9HQq1JipbiGfnmQzgLmK9xdyHPdu/HUJgDCA2QA2cxw3ynGcnGnk0wA+D2BdvOxaAF82KEu4\nEDdkWkiVIXljujpNkxl5jbwKo4DnVs/RIsSL2BXPHbmzl8fLB5IHQznUm1FdvSNCUh19YSEhw/q2\nYbzTM5absmNETHiTSlnl/JXvKcruVNQJINHG+rbhRJCCV9tHUgLDa71P2lt6aCRtPcpMKXK7Zq/h\nLkUuzvVtw1h/MLYRfmdPqvwy6vZXtw4hMJq+73lBwup4EPdlH/bhlYNjm+43xM/T+k6Y+Z644fdW\nLOTtmSTHcU8zxmYo3ncAAGPsDAC/AnA2Yh6hcrd0EMAkABMVx7WOKY8bUls7ASUlqc9W7MLvd88q\nRTvJVq8H3u1JuvtecyiC+Z/R32aQC9QyLNndl5BhNCri9a5uAMDrXaMon1BqSl49vYz0zeQcLa5f\nuz+ROzIqAQdHtAdCvbpGoyLOW7cnqY6rtvSgbWgsIbKaq7f2YM3/zEySVT5fWVpZp/K/ss6WUBS3\n7RpI0ln9fv4bAZx/3OTEsYiYvh5lEAS5XbPXUJmLU9mGoHFMr/0BRUR4o3b/8XZ3omz/qIh+xWet\n8fMkCSnfCa1j6vrd8HuTKfQx0dGFO4yxCwDcAODrHMcFGGODAJQ9WoNYOi75+IjGMXVZQ/r6Mosd\naQa/vwaBQDB9wQIjW714QcLzTcmX5rnGfpw7rTRviwq0ZFizdzAhw9LGYCKtT3M/j54P0surp9eX\n60t09QVg+RytPgrxInYc0g8bp0arroc/GMCQKkOzbCD1OBiKYndLX4qsmbKjc8TwfdewgP97v8dy\nPVqYvYaZYNS+XrvLP+jVPQcAVnF9UF/51Y39gJR8zdT1u+H3JlNIY6KeMXdsCwhj7BLEPMg5HMft\njR9+C8BZjLEKxtgkACcC2A1gM4Bz42XOAbDJoCzhMqwGAM+XDM39fMrUqkw/nzwVqyWvnl6LdvTq\n6pvJOVos3NGLqIUZNHVdvCDhaY04qmb46aYuw8wkVlDroKXToImmzPSF2WtohrkzqvDonKl4dM5U\nzJ1RZdi+VrvLmoJJHqcW+4eERO5QGWU+Ub363fB7KyYcMZKMMR+AexHz/lbGV7gu4jiuM358E4AN\nAG7gOC4M4BYAFzLGNgM4HcDfDcoSLsIN+wPTyaBM9muEUl6jOt9TPNdSsu7AENZZPEerj0K8mHjW\nZwVldpGljYMYzTBWuE320RHMXsN0bO0KY1qlzzB/p1G7z2Z4g5Kufjf83ooNjySNr04LBII5U7iQ\nphaskI1egiil5ExUYiW/YaboyVBXV4VD3SHctL1X8/OpFV7cOntKUpJiWV69OuVsHf6KVL0EUYIH\ngFd13OgcZZsyfFRMGFVRlBAIC+gJC+gYESBBwgSfBzVVFRgZ4TGx1IPach8+UVuGUp8X9eVe/Hpr\nD4b4KA7oPMOsKwPOn1GJlqGxz6tLPfDAA1GSMLHUi9P85RgUxtrv50XIY4kHsdWw8ova8vgaAA/g\nhQdTymP35hKQ0Kuu3IveuNWOCBLu/3AwqX/95R5c9anJmDxpAvZ1hRCVJDy9bwjzT5iIUk+sT0VR\nghTvpq4RAY83BtHLj+kwfYIPi2bV4ej4lhv5GgqihM7haKJ9UZKwsX0Eb/dGUvrmq0dU4sKZ1fB6\nPZg+IaZX+7CAiCjhhrd6NLOXyN8jud2RiIgfbezSLCvr+quTJmOQFzGl3ItpE2Lnyd8TAJrfFaU8\neuTj96akkMZEv79Gs2MomACRU9wQ4kxPBv+kMlSESwzjlOoNKkZ6HT2x1LKMVs4pK/HilKkVhmX0\nBqeljUF8NJA6+Cvp5YEoSvC7kzNbcJHpBvaZ8f9LG4MpA31gVELzYBQLPlWNEyslLG0M4uCQgJZg\nVDOk31KVgQRixqOhM4xj4+WV13CGov95QcKTe7U9vd19PA6vKknS78jqkrTZS5TfozKfJ1FWK31Y\nYFTS1cvM98Tp31uxQb1JjGvcYMTzBS8Y54hUopcv0kwb12zpxuLT6zMylOmmCy8XxJTtDWo509WR\nTq9MArxb+R75vB5Mi0cv0puqzbT/CfsZH6MDQRApOSKVyY5lplb4EtPBmYzPeuHazJLOQHk9nrTp\nxrLNYpLrGyf5RuKuU6cUVXLiYoWMJEGMA7RW8L7bwyeSHdvdhpkUYFrTsmoDpS4jiOnTjbl9dkC+\nkXh631DCuBdDcuJihbKAEMQ4IB/bAvTCtWmRaezRJbv7Cnp7A0XCKTzISBKEi7FjEE33jC7EZ7gX\nRNWGloenJ38msUd5QcKqPQOaZdVbLLTkcwNWbiQId+DeOQmCGOdkuwhGRn5GxwsS/vRuH/735FpI\n0lg6q9++2YO7z8iuDSsJpc1My2qV8XmA+77ysZRnd0o9tfrMrn7MFr0biXkzqx2SiDADeZIE4VLs\nyuIgP6PbciiMvcEo3ugYwd27Ymm6thwKp2SdsIrVDexmvCmtMj6vBzMmlSWlHFP++byejDNm5AOK\nhNo5pTcAACAASURBVFOYkCdJEA6jtYDF6iIYM23I9a1qHUY/L+KJPUFs7Ahn3YaV1aRG3pTS88vE\n49LqM/m1Vjv5xMzWFsKdkJEkCAfRmwpMt83BKsr65Li0srHMtg3latJ0gQTMTMsalVlw2ERTdSs9\nNDv7MVPMbG0h3AlNtxKEg2hNBVpdBKMkXW5BJeog7tmutky3YtXMtGz6Mtoel5aOerkfnVjEI99I\n6P2V5DFUHGEN8iQJwiH0plStLIJR1yd7pUrMZrrI1tNSGnytaU0z07K8IJn2uJReq5aOerkfnfIm\niWQyDV+Yb8iTJMYto1FnnwNpTQ9mk8VBzyu1kuki00wRSoO/vm0YCxpSPcp03pQgAb95swfTKn1p\nPS6l15ovHc3glq0mbsfMPlm3QJ4kMS7hBQk/WncAd86qdWwhh9aU6vnHVGUUUk3tlcoLQWTvzUyW\nitIsw9HJBj/mwQmWPTYrIe2UZbX6TIhnJwGSQ+3J5OKSu2WrSSGQbfjCfEJGkhiXrGgOYVcg7NiP\nVG9KVRmqLNP6WkJRLNndh29PL0t4b1ayVFhF75nnBgurSa2s5lV7rVsOhTX3ec7IIBtLNhTSwO8k\ndq/czjU03UqMO5wODWZ3YlwtI7Vm72BSPemmOrPJMaj3zLPVwh5AK5Fo1F5rtvs87cDp71QhUWhR\nh8iTJMYddm+vsEq2WSrUaBmp5n4+L3qlex5oxps0s3dSZjQqZu215gKnv1OFgpVr7RbIkyTGFdls\nr7ALO706u71Sq8gG/9E5UzF3RlXK52a8SSuRaLQCnJttJ1e44TtVKBRi1CHyJIlxRabbK/JBJkvi\n9bxSOT9hrm/OZYPPC1JGCYStJEg2CnAOOOdNuvk75SayTYbtFGQkCdtw+74nqz/SdPpY0deorOxx\nKPc4aslhJv+ijFZ+wnQyKNswq1uIF1Fd5s14Cll9XkQc875KvZ6k85QBzle3DGFlS3Kg81Ydw2RV\nJ+U5eu+Vx818p7L9bYR4EWU+j6U6nP49pmwBsvkxQ76g6VbCFgph35NyavDROVOxeu4xide3zZ6S\nEmP06oaAYfQYs/oaleUFCVdtDmB5UxAf9kewrCmYVFY+N8SLlvtXWTadDMo2ekcEU22FeBEXrO9E\niBdTppDVex31ppCV502r9OEv7/Xjrnf78Of3+jCtMnnFrc/rweFVsXJbDLxWLb2t9J+6r4z6Tv5O\nPXi2H4/OmYoHzvLj2JoSPHi2P/Gdyva3IffzlZu7TNfh9O9xrP2xLUe5XDyWS8hIErbglkwLRqh/\npOqMEsof6bKmILiBKJY3BTXrsqKvUdmljYNoHIzi2ZbYZ6tbh5LKyucu2tFrqX9Ho8lGwUgG+bOF\n8TYufe2QqbYW7uhFWAAW7ehNOp7pAL2iOYSPBiJoHIyCG0h9TiXvbQ1HRVT6PAnDpPxT3+xk0n/q\nvjLqO5/Xg2mVPizeOZZVZW8wiobOsGFmEivI/dw0KFhaLezk71Fuf8nuPkfatxMykkTWFNvyd16Q\nsLo1NpW3qnUoRR8r+hqV5QUJ/9kXa2cw/khrICIlyob4sZWcO3t5U+3JLNndlxgk08kgf7Yr3sZI\nPJqbUVshXkyU39nLJyVuzmSA5gUpJc7qhrZkr1De2/rHt/uSDJGeR6LUTe6/9W3Gi5nUfaW8Bnr9\nIeu7rCmYUpYXJKxvGzY83whlPwPA+oOp38d0OmSypSiTz7TaV29FKkTISBJZU2j7ntKxrCmYMFYD\nESnFm8x0T5+67NLGQYzqRMZrCUWxcEdv4tyoNHY8Xf/ygoS1+wYBxAbJ5XuCujIo5YuqxjKjthbu\n6E2Uj0pj3mSmA/SK5hBaVYtfWofGPCelETV7w6Cl2/4hQXd2QH2O+hpo9YdS39WtQylllzUFEzFk\nM/ltKPsZAPYPi5ZWC1tt08y0fLprqmxf3opUyJCRJLKi2Ja/K71IGaU3aUVfo7JKL1IPpQdhpj2Z\nFc0hNPfHzm0JRbG6NXlhiVIGrT2H6dpSezfAmDeZyQCt5UXKyN6k0oiauWEw0k1rdkDvHLWe6v5Q\n6ivfWCVkbx/BqpaQ4flGaPUzYOxNZvt7NDMtb3RNi208AMhIEllSiPuejFB6kTJKb9KKvkZljbxI\nGbVnl649QHuQUqfEks83kx1Eqy21dyPLunB7T0YDpJYXKdM6JGBZU1DXiGpNH8t16ummNTugd46R\nd53uJqM1FE1Mo2udnw6tfgaMvclsfo9mp+WNrmmxjQcAGUkiC9Itfy+0u0dekPBsi7Z390zLEEK8\naFpfo75Zd3AYK/Yae5Hp0OtfM4ZPlmHdAXMyKNsK8WJiulPNu70RywMkL0hYlyaDxzMt+kZUbbTk\nlazpsoI80zJk+nqpkfvDbF/rnW+EUT8DwNr9Ic1n5dn8Ho1mAczMEBTbeCBD+ySJjCnUfU96CKKE\nyhIPghq375UlHgDGuQ7Ve/r0yoZ4AVc09Gh+5gWw+Iw6hKOAKEroGRVRV+5NWnnrr4htjVD3r9Eg\npc70IYoSZC3bh6NY/F4f+iKp5/nLPVg0qy7RVpkXuPWUupRyEVHC3Tv70cun9p3RRnGfB7h1Vh06\nh6Po0XCtJ5Z58Pf3BxEySGsm1y9PB67cF8Jts6dAFCV0hQVERAl/2z2QlAGlpswLQZQSF019vZRZ\nRIBYJhEBwO3v9OF/T66FIBr39Y0n1+Lmd/qS2pw+wYdFs+pS9n9qUeYF/vi5Ws0+AQA2uTSljmx+\nj3rTpPNmVidea32mvKZa7ctBLdK172Y8klSY1j1TAoFgzhT2+2sQCOgvCihUxotegiihfTg1Ua9M\nNpkylPBREe8ZeAmfqStDWYn1SR6l/MrBSUZPfjv0zlXfCaKEf34wmBI4AAD+6/AKXHxcDXxeD+rL\nvfjV5m60hKKYUV2Ce86oR3VZrA+XNgbxSGPq9/eHx9foBh7QQq7nh8fX4KKPVxvq+2rbMJY0pXpb\nWm0Czv/GjPoIgOn+U+O0Xlbw+2s0Lzx5kgQRRy96jd2UlXhxytQK2+tVyq8VccfMeXa0rYeRAdJD\nkKAb7u71zjCu/NRkVJd5sbQxeQXvjzd24fEvTQMAS2Hv9PJBqp/JzZtZrasvL0h4uU37WaUbw6+l\nezQAHUfKjbrkAjKSBEHkHCMDZIR6Ck/2kOWwdM+0hHDBzJqU6cDuURHLm4K45Pga01OQRvkgrWT5\n0Aq1J0/TmplqzTdG07TytLxuxCSX6ZILyEgSRAGRiTfmBjJNSKz2UP2TylA6Esbb3TEv+bX2MEQR\nmgtoVrUO4aLjakx5yVqeotLDtJLeSS3z0sZgIvjBvJnVSQbHDdczXzMohQqtbiWIAsHpeJyZoow6\no46ikwlqr+4ZnRXJels9zNSpt7JTJpNtFevbhrGgwVxMWMI9kJEkiALBjnicdg/IZupb0RxKRJ2R\n9z1milbS5ZoyLx48248Hz/bjiAm+pM82dqQPG6g04jLKgAt2bavYPyTgowHjmLBkMN0HGUmCKADs\niI9rt+dipj4tA/Tk3lBSrFcraCVdbh8W0NAZRkNnGG2qFafK0HZq5MwgjzcOJoy4TEsoimVNwZTM\nMUbB1NXoBRvYoBkTViTPUgM39AUZSYIoAOyIj2t3Zggz9Sm9SBleTM0cYgZe0E+6bBQcQcvj4wUJ\nP97YhQ/7I/iPTmCHJ5pDGIlKGad30gs20BqKYpEqJuyS3X2OZ+5wG265aaCntQThcqwuHElXh9Vz\nM62PFyRdwyXHepX3MppBmXRZjdVVmMuaguiOb9SXfVrlZn/l6tnvHz/RlHzKRTjpIvioo+k83zwA\nURjLADMetlakI9PFXnZDniRBuBw74mHananFTH0+D3D6tErN86NSLNycFXxeT0oOUPnv6ImlmDGx\n1JTHpw5iLwdYkqdtp1X6klbPmk0PpfR6lNO0c2dUaeqvZN9AJBF6r9BjndqBHY8X7IKMJEG4GDvi\nYdqdmcFsfUaBAGLnmDNAdqMVxF7m1fYRLG/STy2mh3qqVN5WMa3SZ9gHejhtGJzGTen3aLqVIFyM\nHfFxjTzRTKaxzNYniGOxbuXpSyWtWciQKVqp0JS0hKLoU32ebvrTaOpZff3keLIy/gofXtg/nNI3\n2VyfQseOxwt2Qp4kQbgY2SPJZOEIYH9mBrP18YKE37zZg2mVPkNvKt/ZIZYbeJEy6s+zSaytvn5H\nTyzFKVMrEn/Tq0pc0zduwW3ptsiTJIgixu5MLWbrU04/XvTxaldki4nFVNU28JU+YEQAJvgArbjl\nenFKs/V65P7UCkgvfz6eSHcT5oQ3SUaSIIoYu0OOmQ1kbjYYeD7xeYDbT61POR4RJSzc3ou2YQG1\n5T7cE1/hqnW+mmynsuX+tBKQvphxY/o957+5BEEUFVaCgecTPQO/tDGYCELQFl/hakZeN3o9hY4b\n48i6SxqCIAoaty26SEc28rrR6yHshxbuEARhG25bdJGObOTNdlEVURiQkSQIuCNGZKGT7UrafF8D\nu1f+EsUJTbcS455MEwITyWQz/ejENaDpUsIMZCSJcY9bYkTmhAgPlJYl3uYyya8gIeNFF+proI6D\narSR34o+yvLqRSKy56hs15daRdaoZdbyWNWyOIUd3xctfTO9Zk6Q1+lWxtipjLHX4q8/zhh7gzG2\niTF2P2PMGz9+E2PsLcZYA2NsttWyBGEFN8WItJ0ID99fbgciEQC5zaqQTd3qayCnsJLzOerVa7XN\ndHVd3RBIJEXOVV+p6x17LyaVUcriFHb0wWhU1NE3+2uWL/JmJBlj1wF4CEBF/NBiAL/nOO4sAB4A\n32KMfRbAFwCcCuBCAPdlUJYgTOOmGJF243lpLTz7muF5eS0A+1NlKQeubOpWX4Mbt/ck6jKq12qb\n6eriBqKJpMi5Slulrld+v2R3n64sTmFHHyzZ3aeprx3XLF/k05NsBjBX8f5zADbGX68F8GUAnwfw\nEsdxEsdx+wGUMMb8FssShCm0stwXizcp8Ty827cCALzbtoAP8xl5zHrllHf42XjjWlsw3uuNeb7r\n24axQadeq20qy69vG07ILX+mTAxt1G42KGVQJ15es3cwIZNSlg05/j4aXV87knyv3TeYqCM10bT5\na+bk7zJvzyQ5jnuaMTZDccjDcZysdRDAJAATAfQoysjHrZQNGMlRWzsBJSW5eNIQw+8vsmdacYpR\nrwfe7dFc/r/mUATzP6O/oKMQGHnqKXja2wAAnvY2ND/7PFqqTwNgXsfRqIgfrTuAR845EmW+5Pvp\nB96NeXtrDkUgSUjyBK30n9Y1kFEma07Ue1jsu6g8z0ybyvL7hwQ8czCMLZ0jeOScI/Horr6ktjTb\nteH7oJShNRTFLTv7E++b+/lEXyrbb83h9zHd9c30mirraO7nE3Xctmsg42vm5O/SyYU7ouJ1DYB+\nAIPx1+rjVsoa0tennwg1W/z+GgQCwZzV7xTFqJdRlvunP+rDudNKHV80kTERHmVvbEo6NPntLSg9\ncxYivthP/rnG/rQ6Lm0MYlcgjPu2diYtaOIFCc83xX5qqxv7ASn5Dt9M3XI9K7k+wzLqei87qRaB\nQCjRvpk2lfLKLHm/D6MicE9DBzZ2jKScY7Zus2jJsKMzObD56sZ+iKIINattaF8LM9dXxmofaOub\n3M9Wr5kd18EIPUfAyX2S7zDG5sRfnwNgE4DNAL7GGPMyxo4C4OU4rttiWYJIi5zl/tE5U1P+bps9\npaCX/3teWgvxwIGkY0cOHMIFza8n3qd7/mo01aV8htgaiqJ1KDkiuOnN+IrExI/OmYpPTCo1LN8S\nimLJ7j7LAQC0yo/GbdHq1qEkz02vXTuTVMuozWFrKIoDw6lGstWG9tWYvb4y2ST5llEnmrZ6zZxa\nM+CkkbwWwCLG2BYAZQD+w3HcDsQM4BYATwO4PIOyBJEWoyz3BR0tJRKBd8smzY/++8B2lApjA4/R\nhnm9BU1azxC1MLMZXxmxprbMi8bBSNp6n2nsx7o0AQDUzy71AgYAqWmx9MgmuEA6GcywzubgBkbX\n144k32b11arTbUEePJJU+IsUrBAIBHOmcDFOSwKkV0EhCEB3VyL1kigChxRJfqN1fsA39kx++gRf\nyg0BL0j4xaZA0p38jOoS3H+WHyuaQ3ikMbXP5s6ownkzqpKOadWtx2PcIJY06XsJ35kxAd+aUY1J\nkyegr3cIXp1668u9+PXWHtx9RiwogSBKaB8WwAsSFu7oRbtGHqypFV5ccdIkeD1Ab3jMk5tS7sVh\nE0oSbVnRR4ksgzLbiLLtP86qw9T6anxwcAC9oyLqKrzweZLbmVrhw8dsunkzur4+DzT7SMZMH8j6\nAkBdXRW6u0MpiaaVdajrVJ6fqQyZ4PfXaFZKwQQIopjw+YBph8PnrwFKg/ACONxiFXpTXcubgnhF\nx4vc2hXG/BMnZvS8KJbn0dg7faebx09P9OGI2nIEorxuuSXcID4aiOCJPUFcyiYmvNWljUHdgbcr\nLGLvYDRngSSUMrSpZOgKi9jaNYoFx03BxGhlTtpXky69V7ZZOJRBGmIpwEpx9ETjqXS9892AeyQh\nCMJxjKa6Xm4bxs2nTNHMtQhkHsZNGR5udcsQVrakJh+WB/EFh03UrYcXJKxujcm+qnUYF368BmU+\nj6npv1yntko3hXi5kPos0gk53JqtxUnISBIEkSBdPNNcTHXJngMvSNgRCOuWS2dMlu8Jop+Pfd7P\ni2PeZFwnUZSSpv1k5Om/XNqGdP3q9eTHMFG8WuuQkSQIIoGTU10rmkMpK2WB5OedesZE6UXKKL1J\nWScr0352kq5fS/K0UMxtU5mFAKXKIgjCcYxWzb7dPYpplT4cWV2ia0yUXqSM7E0SRDaQkST+f3t3\nHyRZVd5x/Ns9L/s6bJbdWRAK3QTdY9QIIlEgAbcSXgSsaJkSKUOZKEURQ0xQqmLF+IIUSkwAExRM\nQgRRAgYVsiXlvhiMy4qrMFEiIJ59gV2WWWCH3Z1hZmdnuud25497e/dO973dt6e770vP71NF0dN9\n+t7z3Omdp8+9555HJHFBk0kqotzT+UDAdUyA+3cdankJviS02pc0xVJPO/rZ6ViVJEUkUQWnzIY9\nwUmuot79cU6pzKKQi2mLvNtAovQh6WoT7epLmmKppx39jCNWnZwWkUT15ODM4xZx/65DnH/iIv7g\nRPdWiIeGD/OD4cO85zWLec9vLg2dVNLfk+MfzlgZuv0oszXTVFO01b6kKZZ62tHPOGLVSFJEEuWU\n3euOANvGipyyYgGnrFjA9jF3BZ7H9xc4blH4rFr/yj1zWT0pLdUm2tGXNMVST7uqjMQRq5KkiCQq\naIm0OOt8pqmmaKt9SVMs9bSjn3HFqiQpIokJmtX6w72H+eHw7Ns5OjVSCNp/vX11urZjK/VNW31/\nXNrRzzhjVZIUkcQEzWptpbJIO/Yftq9OTxJptfJFmipn1NOpKiOdilVJUkQS0Wx1DHeGa/uWb2u2\n2oR/kki7FZxSS5Uv0lY5I0ynq4x0IlbNbhWRRAQtkeZfOq66WgS0d/m2ZpZoq54k0u41TvO5XN2+\nOKVy3TXjsrLcXDv6GXesSpIikoiwJdLqLR3XzuXbmlmiLWiSSDtvOeit05fKad6bz1wZmpizstxc\nO/oZd6w63SoiLUvL6bxOSHpCTCdP80pjSpIi0pKsrPAyV0lOiJkolDJx32M3U5IUkZZ080gnyQkx\nBafM5Zv3ZeK+x26W/pPYIpJanZ7QkrQkJ8Tcs32cl6dnz+btxmOcdhpJSncpFpLuQfOy2GdP6Kon\nXkztHGn5txXllohQ1ce7zvFvdcm7uXLrY9Yu+r5rYiaw/NdcFz+o9/pclsSrNj0z91t20nJqWUlS\nukexQM9NN0CxmHRPostinz2hE1qmpum56QYKU4W2Xas8et2z1PAaaN3Xq493So//vdvHGSsGx3fP\njgkmfLUzq+Ot/H+iUJrzcYp6ndm/z+r2BafMhzfsoeCU55RwG+0/riSqJCldI7dpPblnd5L7wfqk\nuxJZFvtcETah5df3rSP37E6evm9d265VVq573vXkwYbXQOu9Xn2803j8C06ZzS/M/vJx4uIebj9n\nkPeuXkKxDPc/ezQ2f7yV5DJRKHH55n1zPk5RrjP7E1lQ+/t2TvDEyBTf2jHe9JelRvuPc7KYkqR0\nh2KB/NBPAcg/tjV1I4NAWeyzJ2xCS79TZOUTjwGw6slH6XNmeGi4tQku/hHrgzvH+GGd2Z51K0NU\nH+/JydDjn+Spvvt2TtQsyzc86fDjFw4fqZay+YWpIyM0f7z37hjn6dEinxnaf+R6ZrPHKWp1jUoi\nu2f7eE17/zbW7Z5s6stSlP3HOVms59prr+34TtJkcrJwbae2vWTJAiYns3t9KUwW4spteJD80M/c\nx+PjlPv64XWm7nuSjmsufY4qjtjetmoh71m9ZNZ/l+3czPFPDQGwbPoQ0/lethyzmh7g1JUL5rSf\ne3dM8KMXpgAYnS4x5p1qHC2U6M/nePOKBcFtq16vOd67nyH/7DNHf/aOf2WUcu6Jizt2zbFa5fdV\ncMrc/MQo4wGnWp8ZL7JvanbsTx0szIp398QMU06ZkcMlKlto9jjVe62i4JT5ylNj7j7Hi7xUp19T\nXpIbnS5x4UmNj2mj/fv3HXWbUSxZsuBzQc9rJCnZ5xshVKR+ZJbFPvsETmhZUGL5Lx+d1e4P9z5O\nnzPDut2HWq4ZGMQ/0qh703/Q8d6xbfbP3vFP8paWymzar69dNeu/288Z5Ji+2X+ug6qljHpfIKqn\ny0Q9TlEXTvCfaq++dvrQ8OSR0b5flFtYouw/7nJgSpKSeblN68ntHZ793N7hVF1nqpbFPjcSFNNv\nje/j/TsfZqxY5t7ttbMyGwm67unn/yNZ76b/wOM9M7ttbu8wpY3fj+Xm/bDths2m/cmLUwxPzj4F\nG1QtJUzU4xRl4YRGX1yeO+SwO+R31uiYNtp/EqsfKUlKthWL5LduCXwp/5OH0zkyy2KfG6kT0zv3\nDNHnzPDAruZGk1GrhGx8fpKJQngVjYd2v0IupG/VprY8zPCYe6qvU6OUZiedNFstJUyj47Th+Uk2\n7Km97aTy3kp/G31xqafeMY2ycEMSqx9pMQHJtnwe56PX1H09dbLY50a8mEoleGnKoeiU+fJTY+w7\n7I50SrkcA/35htUs/Pw38julMiNTDsuWLWZszP1DumphD3nvWlR/vs5N/47DzJuuOXpYHQcOHjj6\n+vIV0JOn6MD1QwcocbR/nbh53386N8oi6Y2qpfzGgjyfeezAkeuVQVYtzPO504+te5xKpTJlCL2+\n15Orn8hWLcxz3e8ey+h0iWKpzJefHAvs08bnJwOPaaOFG5xS/STaqUUWcuVyOm7YjMvIyHjHAh4c\nHGBkpPlTSmmnuLIn6dju3jbOndtq9/+hNQMtVc/oZFyd6rNfwSnzkS0j7JqYYfXSXr569iD9PbmW\n4nJKZfZOOrMSp1+l5NgJi2tLj811X2Eq+6i0O/bYJRw4UDs6nUtfou57rgYHBwLfrJGkiLRVo9Nm\naVxWLa4+d6Lklr90VL0yY+0QtUxVpd3gsn6WFKZj3Xe7KUmKSFtlpQCwXxx9Dpt0csnJS1vfuHSM\nkqSItFVWCgD7xdHnepNOPnb8MR3dt8xdBmcIiIhkS+OZm3NfCFw6K1tf90QkFQpOOXXXFdOs0enc\nfK47j2UrVUCiiONzqJGkiDQlzsWlu0Wjklu9MS1/Fyd/FZBObT+Oz6GSpIg0Jcll2yQ7KlVAOvU5\nietzqCQpIpFFrRAh81unPydxfg6VJEUksrgXl5ZsqLcAeic+J3F+DpUkRSSSJBaXlvSrvjbY6c9J\n3J9DJUkRiSSJxaUl/aqvDXb6cxL351BJUkQailKhQeaf6muD9aqMtONzksTnUPdJikhDWVxqTjqv\n+trgA7smZn1Oqhc4b/VzksTnUElSRBrK4lJz0llB1wZ/tHeK9588cOQG/3YucA7JfA51ulVERJo2\nX65RK0mKiEhT5tM1ap0/ERGRpsyna9RKkiIi0pT5dI1ap1tFRKShTi5UHte+5kJJUkTmnTT9Ec6C\ndlTcCEuG1dtNW5WZRMfLxpg+4C5gNeAAVwAzwNeBMvAkcJW1tmSM+Sxwsff61dbaR40xrw1qG3MY\nIpIhlT/CN5+5UjUxI/KvqnPZmoGm3x92zIO22+q+2i3pkeRFQK+19izgOuDzwM3Ap6y1ZwM54N3G\nmNOAdwBvBy4FbvXeX9M25v6LSMao1Jcr6kitHRU3go550HbTWGUm6Suv24BeY0weOAYoAmcAm73X\n1wPnAxbYZK0tA88ZY3qNMYPAWwPaPlBvh8uXL6a3t6ftgVQMDib/zacTFFf2dGtsrcQ1PVPi4X0v\nA/DwvmmuOuN4+nuSHiu44vx9Tc+U+PCGPdx54UkN4/+3x/fPWlXn+y8VueKU8Jmt1Y5ZviTwmAdt\nt1ympX11QtJJcgL3VOuvgZXAu4BzvGQIMA4sw02g+33vqzyfC2hb18GDwff2tMPg4AAjI+Md235S\nFFf2dGtsrcZ197Zxdo4WANg5WuDWn76YilN6cf++7t42zhMjUw3jLzhlHtw+Ouu5720b5aLj+iKd\nqh4cHOC2n71Yc8wvOXlpzXbX2YPkcrO32cy+WhX2JSXpr1AfAzZaa9cAp+Ben+z3vT4AjAKveI+r\nny8FPCciUkOlvlzNnNJsdVWd6ZlS4DG/d/t4zXafO+SwO4Ur+CSdJA8CY97jA0Af8AtjzFrvuQuB\nLcAjwAXGmLwx5tVA3lr7ckhbEZEa82UZtUaiFixux6o6dz15MPCYP7DrUMg75r6vTkn6dOuXgDuM\nMVtwR5CfBIaA240x/cDTwHestY7XZituYr/Ke/811W3jDkBE0q/RH/xLTl46L2a6ho2mg+JvdVWd\nglNm3Y6xwNcW9ea4+ayV9OXdjTilMiNTDgCrFvaQz9f2JSm5cnl+nWoYGRnvWMC6DpQt3RoX/P0k\nmgAACTNJREFUdG9sc43LKZXZO+mEvn7C4h568tH+EheccvsSarEAff2x/b7u3jbOndtq9/OhNQNt\nvzbrlMpMLVw4q1SWXzPHPA6DgwOBnUn6dKuISMdVllEL+6+ZBNm2G92LBXpuugGKxda3FUHci5L3\n5HOsXtbf8jFPWtKnW0VEMqOdN7rnNq0n9+xOcj9YD3/6J23qYbj5tCh5O2kkKSISQVtvdC8WyA/9\nFID8Y1spxzCabNdoer5RkhQRiSDqrNAocpvWk9s77D7eO8zUuv9qSx+l/ZQkRUQaaOUey5o2vlHk\nkTZbtsR2bVKaoyQpItLAXO+xDJro4x9FVpT27HGvTUrqKEmKiNTRyqzQmoW9i0XyW4PXPMn/5GGN\nJlNIs1tFROqY66zQ6ok+l5y8lP58Huej19S0PfbYJe79hHmNW9JGSVJEpI7KrNBmBU30uWzNABz3\nqtp9DA5AX/ct/tAN9LVFRKTNtJh691CSFBFpMy2m3j2UJEVE2iju5d+ks3RNUkSkjbT8W3dRkhQR\naaO5TvSRdNLpVhERkRBKkiIiIiGUJEVEREIoSYqIiIRQkhQREQmhJCkiIhJCSVJERCSEkqSIiEgI\nJUkREQk135fRU5IUEZFABafMx7e+PK8TpZKkiIgEum/nBE+PFud19RIlSRERqeGviTmfa2EqSYqI\nSA1/Tcz5XAtTSVJERGbxjyIr5utoUklSRERm8Y8iK+braFJJUkREjig4ZTY+Pxn42sbnJ+fdaFKV\nQUVE5IieHHzhbSvqvj6fKEmKiMgRPfkcJy1VaqjQ6VYREZEQSpIiIiIhlCRFRERCKEmKiIiEUJIU\nEREJoSQpIiISQklSREQkhJKkiIhICCVJERGREEqSIiIiIZQkRUREQuTK5fm1oruIiEhUGkmKiIiE\nUJIUEREJoSQpIiISQklSREQkhJKkiIhICCVJERGREEqSIiIiIXqT7kBaGWP6gDuA1cAC4HrgV8DX\ngTLwJHCVtbZkjPkscDEwA1xtrX3UGPPaqG3jjKvCGLMK+F/gPK8vkfqa5riMMX8L/BHQD9wGbI7a\n17TG5X0O78L9HDrAFWT892WMeTvwRWvt2mb61462McZ1KvBl3N/ZNPBBa+1LxpgrgCu9vl5vrX3Q\nGLMSuAdYBOwFPmStnWymbVxx+Z77APBRa+2Z3s+ZiysqjSTDXQbst9aeDbwT+ApwM/Ap77kc8G5j\nzGnAO4C3A5cCt3rvb6ZtrLw/vP8KHJ5DX1MZlzFmLXAW8HteX05qsq+pjAu4COi11p4FXAd8vsm+\npiouY8zfAP8OLJxD/1pqG3Nc/4ybRNYC9wOfMMYcD/wV7mf0AuAGY8wC4DPAPV5ffwFc2UzbmOPC\nGPMW4HLc40oW42qGkmS4bwOf9h7ncL/1vBV3dAKwHjgX+H1gk7W2bK19Dug1xgw22TZuNwL/gvuN\nDbojrguAJ4AHgO8BDzbZ17TGtc3bbx44Big22de0xbUTeK/v507FEtS2k6rjutRa+7j3uBeYAt4G\nPGKtnbbWjgE7gDd7MWyo6mszbTtpVlzGmBXAF4CrfW2yGFdkSpIhrLUT1tpxY8wA8B3gU0DOWltZ\nx28cWIb7h2vM99bK8820jY0x5s+AEWvtRt/TmY8LWAmcDrwP+HPgP4B8F8Q1gXuq9dfA7cAtZPj3\nZa39Lm6ir+hULEFtO6Y6LmvtCwDGmLOAvwS+VKev/ucbxRXUtmP8cRljeoCvAR/39l2RubiaoSRZ\nhzHmJOB/gG9aa+8B/Nc0BoBR4BXvcfXzzbSN04eB84wxPwJOBb4BrAroU9bi2g9stNYWrLUW95u7\n/x9aVuP6GG5ca4BTcK9P9gf0KWtxVXTq31RQ21gZY96Pe8bmYmvtCOF99T/fKK6gtnF5K/A64KvA\nt4A3GGP+iezHVZeSZAhjzHHAJuAT1to7vKd/4V37ArgQ2AI8AlxgjMkbY16NO3p5ucm2sbHWnmOt\nfYd3reRx4IPA+qzHBfwYeKcxJmeMOQFYAjzUBXEd5Og37ANAX5N9TWtcFZ2KJahtbIwxl+GOINda\na5/xnn4UONsYs9AYswz4bdxJRY/gXnv297WZtrGw1j5qrX2j97fjUuBX1tqrm+xr6uJqRLNbw30S\nWA582hhTuTb518Atxph+4GngO9ZaxxizBdiK+6XjKq/tNcDtEdsmrZm+pjIub4bcObj/CCt9eLaJ\nvqYyLtzTdHd4/ejH/VwONdHXtMZV0anPXk3buALyTkveAjwH3G+MAdhsrf2sMeYW3ASQB/7OWjtl\njLkeuMu4sz5fBj5grT0UtW1ccYWx1r7YjXFVqFSWiIhICJ1uFRERCaEkKSIiEkJJUkREJISSpIiI\nSAglSRERkRBKkiIiIiGUJEVEREIoSYp0EWPMoDFmvfd4lTFmOOk+iWSZVtwR6SLeGqEXej+ejlsz\nVETmSElSJAOMMV8ETsCtdvIG4L+Bb+JWp3k9cLu19nNeu33W2ptwk+SQ9/7lwN8DvwOsAL5rrf2k\n99oZuOXTFuMuLv0la+1t3murcAtYrwFewi3dtcta+4/1tinSLXS6VSQbTgMW4BYPfhPu2pYXAucD\nZ+OuK1xp93Pv8enAkDEmh1u14dte8eY3Aucbt3gugMVdiPs07z2fB/Detw5Yb619M/AXwEciblOk\nKyhJimTDW4BrrLUF3ALgDnCdtbZSEuqA9/9TcSu7g1vaaAg4D7fY7Y3GmMe951YCPV67i3Erpvwf\nbjWGyjbPxa3L+DUAa+12YBL3FG6jbYp0BZ1uFUk5Y8xrgHFr7R7vqVOAX1prD3k/nwb83Gv3irV2\n1CsX5lhr93mju68GnQo1xrwLd3T4x141h/cBV/q2O+Rra4Bha+0r9bYp0k00khRJv9OYPQGnMkKs\nfr3mVKv3+HncQtuLAYwxS7yEd+S9XoI8Ebje974R4PVercZFuOWfomxTpGsoSYqkX5Qk+XPCk+S3\nvMe/9E6NPgS8ynvtm8Ba71TrjbjJr7Kv/8Q9rWuB+4EJ772NtinSNVRPUkQCGWMGrLXj3uNTgW8A\nZ1trx5LtmUh8dE1SRMJcaYy5HDgM7AcuU4KU+UYjSRERkRC6JikiIhJCSVJERCSEkqSIiEgIJUkR\nEZEQSpIiIiIhlCRFRERCKEmKiIiE+H8SfSBEkcoOlwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10814e650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "outlierIndex = train_error >= np.percentile(train_error, 95)\n",
    "inlierIndex = train_error < np.percentile(train_error, 95)\n",
    "\n",
    "# 得到train_error最大的index值，即极端离群值\n",
    "most_severe = train_error[outlierIndex].idxmax() \n",
    "\n",
    "fig = plt.figure(figsize=(7,7))\n",
    "indexes = [inlierIndex, outlierIndex, most_severe]\n",
    "color = ['#2d9ed8','#EE5150','#a290c4']\n",
    "label = ['normal points', 'outliers', 'extreme outliers']\n",
    "for i,c,l in zip(indexes,color,label):\n",
    "    plt.scatter(train['mileage'][i], \n",
    "                train_y[i], \n",
    "                c=c,\n",
    "                marker='^',\n",
    "                label=l)\n",
    "plt.legend(loc = 'upper right',\n",
    "           frameon=True,\n",
    "           edgecolor='k',\n",
    "           framealpha=1,\n",
    "           fontsize = 12)\n",
    "plt.xlabel('$mileage$')\n",
    "plt.ylabel('$price$')\n",
    "plt.grid('on')\n",
    "sns.set_style('dark')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们来看看离群值的数量有多少？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False    396\n",
       "True      21\n",
       "Name: price, dtype: int64"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "outlierIndex.value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上图结果也符合我们的经验理解，二手车的行驶公里数越高，它卖出去的价格就应该越低，所以对于处在右上和左下区域的点可能是一些离群值（对于同一款车而言）。比如左下区域的点，一些行驶里程数低，价格也比较低的车辆，有可能该车辆是事故车辆或者有损坏，而右上区域的离群值有可能是真实的离群值，相对来讲不容易有合理的解释，可能是输入失误或者胖手指输入造成。\n",
    "\n",
    "本案例中的数据只有400多条，如果数据再多一些，则检测的结果会更加可靠。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. 标准化对离群值检测的影响"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通常情况下，为了避免不同尺度的影响。我们在进行线性回归模型拟合之前，需要对各个特征进行标准化。常见的标准化有z-score标准化、0-1标准化等，这里我们选择z-score标准化来观察标准化对离群值检测的影响。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'LR_score' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-12-f0395cb1470f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      8\u001b[0m ' + '.join([format(LR_zscore.intercept_, '0.2f')] \n\u001b[1;32m      9\u001b[0m            + map(lambda (f,c): \"(%0.2f %s)\" % (c, f), \n\u001b[0;32m---> 10\u001b[0;31m                zip(dv.feature_names_, LR_score.coef_)))\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'LR_score' is not defined"
     ]
    }
   ],
   "source": [
    "# 利用 preprocessing.scale函数将特征标准化\n",
    "columns = train_x_dummies.columns\n",
    "train_x_zscore = pd.DataFrame(preprocessing.scale(train_x_dummies),columns = columns)\n",
    "#train_y_zscore = pd.DataFrame(preprocessing.scale(pd.DataFrame(train_y,columns=['price'])),columns = ['price'])\n",
    "\n",
    "# 线性模型拟合\n",
    "LR_zscore = LinearRegression().fit(train_x_zscore.values,train_y)\n",
    "' + '.join([format(LR_zscore.intercept_, '0.2f')] \n",
    "           + map(lambda (f,c): \"(%0.2f %s)\" % (c, f), \n",
    "               zip(dv.feature_names_, LR_zscore.coef_)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_y_zscore = LR_zscore.predict(train_x_zscore)\n",
    "train_error_zscore = abs(pred_y_zscore - train_y) # 计算绝对误差\n",
    "np.percentile(train_error_zscore,[75,90,95,99]) # 计算绝对误差数据的百分位数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outlierIndex_zscore = train_error_zscore >= np.percentile(train_error_zscore, 95)\n",
    "inlierIndex_zscore = train_error_zscore < np.percentile(train_error_zscore, 95)\n",
    "diff = (outlierIndex_zscore != outlierIndex)  # diff 用于存储标准化前后的离群值检测结果不同的index\n",
    "diff.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# 画出标准化前后的检测差异点\n",
    "fig = plt.figure(figsize=(7,7))\n",
    "\n",
    "# rep_inlierIndex为标准化前后都为正常值的index\n",
    "rep_inlierIndex = (inlierIndex == inlierIndex_zscore)\n",
    "\n",
    "indexes = [rep_inlierIndex, outlierIndex, outlierIndex_zscore]\n",
    "color = ['#2d9ed8','#EE5150','#a290c4']\n",
    "markers = ['^','<','>']\n",
    "label = ['inliers', 'outliers before z-score', 'outliers after z-score']\n",
    "for i,c,m,l in zip(indexes,color,markers,label):\n",
    "    plt.scatter(train['mileage'][i], \n",
    "                train_y[i], \n",
    "                c=c,\n",
    "                marker=m,\n",
    "                label=l)\n",
    "plt.xlabel('$mileage$')\n",
    "plt.ylabel('$price$')\n",
    "plt.grid('on')\n",
    "plt.legend(loc = 'upper right',\n",
    "           frameon=True,\n",
    "           edgecolor='k',\n",
    "           framealpha=1,\n",
    "           fontsize = 12)\n",
    "sns.set_style('dark')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从结果可以看到，绝大多数样本的检测结果一致。有两个样本存在差别，其中一个样本在标准化之前会被检测为离群值，另外一个样本在标准化之后会被检测为离群值。虽然在本例中，标准化前后的检测效果差异不是很大，我们仍然建议在线性建模之前对特征进行标准化。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. 测试集的验证"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们先以$mileage$为横坐标，$price$为纵坐标画出训练集和测试集的所有样本点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = pd.read_csv('accord_sedan_testing.csv')\n",
    "\n",
    "datasets = [train,test]\n",
    "color = ['#2d9ed8','#EE5150']\n",
    "label = ['training set', 'testing set']\n",
    "fig = plt.figure(figsize=(7,7))\n",
    "for i,c,l in zip(range(len(datasets)),color,label):\n",
    "    plt.scatter(datasets[i]['mileage'], \n",
    "                datasets[i]['price'], \n",
    "                c=c,\n",
    "                marker='^',\n",
    "                label=l)\n",
    "plt.xlabel('$mileage$')\n",
    "plt.ylabel('$price$')\n",
    "plt.grid('on')\n",
    "plt.legend(loc = 'upper right',\n",
    "           frameon=True,\n",
    "           edgecolor='k',\n",
    "           framealpha=1,\n",
    "           fontsize = 12)\n",
    "sns.set_style('dark')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们来看看利用在训练集上训练得到的模型在测试集上的泛化效果："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pred_y_test = LR.predict(dv.transform(test.T.to_dict().values()))\n",
    "test_error = abs(pred_y_test - test['price'])\n",
    "\n",
    "# 使用分布图观察测试集误差\n",
    "fig = plt.figure(figsize=(7,7))\n",
    "sns.distplot(test_error,kde=False)\n",
    "plt.xlabel('$test\\_error$')\n",
    "plt.ylabel('$count$')\n",
    "plt.grid('on')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 找出极端离群值\n",
    "most_severe_test = test_error.idxmax()\n",
    "test.iloc[most_severe_test]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从分布图中可以看到，我们的模型对测试集上其中一个样本的预测表现非常差。该样本是一个极端的离群样本。该车是一个6缸高配版的车，并且其已行驶英里数只有~60,000英里左右，但是其卖出的价格才\\$2612。\n",
    "\n",
    "根据经验，我们猜测这个离群样本出现的两种可能：\n",
    "\n",
    "1. 在网站里填写时出错；\n",
    "2. 该车辆有车体的损伤或者有汽车所有权问题（偷来的或者劫来的）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. 在测试集上使用LOF进行离群值检测（待在样本点上画出代表LOF的圈20170512 18:30）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.1 原测试集数据可视化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = pd.read_csv('accord_sedan_testing.csv')\n",
    "fig = plt.figure(figsize=(7,7))\n",
    "plt.scatter(test['mileage'], \n",
    "            test['price'], \n",
    "            c='#EE5150',\n",
    "            marker='^',\n",
    "            label='testing set')\n",
    "plt.xlabel('$mileage$')\n",
    "plt.ylabel('$price$')\n",
    "plt.grid('on')\n",
    "plt.legend(loc = 'upper right',\n",
    "           frameon=True,\n",
    "           edgecolor='k',\n",
    "           framealpha=1,\n",
    "           fontsize = 12)\n",
    "sns.set_style('dark')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.2 计算每个样本点的LOF因子值"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 可达密度计算公式\n",
    "\n",
    "$$lrd_k(x) = \\left( \\frac{1}{k} \\sum_{y \\in N_k(x)} rd_k(x, y)  \\right)^{-1}$$\n",
    "其中 $rd_k(x, y)$ 样本点$x$到样本点$y$的第$k$可达距离"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LOF因子计算公式\n",
    "$$lof_{k}(x) = \\frac{1}{k}  \\sum_{y \\in N_k(x)}  \\frac{lrd_k(y)}{lrd_k(x)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LOF算法的一般流程可以描述为：\n",
    "\n",
    "1. 初始化$k$，用于后续计算第$k$距离；\n",
    "\n",
    "2. 计算每个样本点与其他点的距离，并对其排序；\n",
    "\n",
    "3. 计算每个样本点的第$k$距离，第$k$领域；\n",
    "\n",
    "4. 计算每个样本点的可达密度以及LOF值；\n",
    "\n",
    "5. 对所有样本点的LOF值进行排序，与1作比较，越大于1，越可能是离群值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_2d = test[['mileage','price']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 方法一"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.neighbors import NearestNeighbors\n",
    "import numpy as np\n",
    "neigh = NearestNeighbors(5) # 默认为欧式距离\n",
    "model = neigh.fit(test_2d)\n",
    "\n",
    "data = test_2d\n",
    "# dist为每个样本点与第k距离邻域内的点的距离（包括自身）,neighbor为第k距离邻域点的编号（包括自身）\n",
    "dist, neighbor=neigh.kneighbors(test_2d,n_neighbors=6)\n",
    "\n",
    "k_distance_p = np.max(dist,axis=1)\n",
    "\n",
    "nums = data.shape[0]\n",
    "lrdk_p = []\n",
    "lof = []\n",
    "for p_index in xrange(nums):\n",
    "    rdk_po = []\n",
    "    neighbor_p = neighbor[p_index][neighbor[p_index]!=p_index]\n",
    "    for o_index in neighbor_p:\n",
    "        rdk_po.append(max(k_distance_p[o_index],int(dist[p_index][neighbor[p_index]==o_index])))\n",
    "    lrdk_p.append(float(len(neighbor_p))/sum(rdk_po))\n",
    "\n",
    "for p_index in xrange(nums):\n",
    "    lrdk_o=[]\n",
    "    neighbor_p = neighbor[p_index][neighbor[p_index]!=p_index]\n",
    "    for o_index in neighbor_p:\n",
    "        lrdk_o.append(lrdk_p[o_index])\n",
    "    lof.append(float(sum(lrdk_o))/(len(neighbor_p)*(lrdk_p[p_index])))    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 方法二（需要检查一下，有一点问题）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 导入sklearn用于计算最近邻相关数据的方法 NearestNeighbors\n",
    "from sklearn.neighbors import NearestNeighbors\n",
    "import numpy as np\n",
    "\n",
    "data = test_2d.values\n",
    "\n",
    "# 定义函数计算第k距离及其第k距离邻域内的点\n",
    "def k_Distance(data, k):\n",
    "    neigh = NearestNeighbors(k)\n",
    "    model = neigh.fit(data)\n",
    "    \n",
    "    nums = data.shape[0]\n",
    "\n",
    "    k_distance = []\n",
    "    neighbor_info = []\n",
    "    \n",
    "    for index in xrange(nums):\n",
    "        \n",
    "        # dist 包含第k距离邻域及其内的点与目标样本点的距离信息\n",
    "        dist = neigh.kneighbors([data[index]], n_neighbors = k + 1)\n",
    "        # 最大的dist\n",
    "        k_distance.append(dist[0][-1][-1])\n",
    "        \n",
    "        # neighbor为k个邻居的索引，dists为相应的距离\n",
    "        dists, neighbor = neigh.radius_neighbors([data[index]], radius = k_distance[index]+1)\n",
    "        \n",
    "        # 排除自身\n",
    "        mask = neighbor[0] != index     \n",
    "        neighbor_info.append([neighbor[0][mask], dists[0][mask]])\n",
    "    return k_distance, neighbor_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 定义函数计算局部可达密度\n",
    "def reach_density(data, k_distance, neighbor_info):\n",
    "    \n",
    "    nums = data.shape[0]\n",
    "    density_list = []\n",
    "    \n",
    "    for index in xrange(nums):\n",
    "        \n",
    "        neighbors, dists = neighbor_info[index]\n",
    "        nums_neigh = len(neighbors)\n",
    "    \n",
    "        sum_dist = 0\n",
    "        \n",
    "        for item in xrange(nums_neigh):\n",
    "            \n",
    "            k_dist_o = k_distance[neighbors[item]] \n",
    "    \n",
    "            direct_dist = dists[item]\n",
    "        \n",
    "            # 计算p点到o点的第k可达距离\n",
    "            reach_dist_p_o = max(k_dist_o, direct_dist)\n",
    "        \n",
    "            sum_dist += reach_dist_p_o\n",
    "           \n",
    "        density_list.append(nums_neigh/sum_dist)\n",
    "    \n",
    "    return density_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 定义函数计算LOF因子\n",
    "def cal_lof(index, data, neighbor_info, lrd_list):\n",
    "    \n",
    "    point_p = data[index]\n",
    "    neighbors, _ = neighbor_info[index]\n",
    "    \n",
    "    nums_neigh = len(neighbors)\n",
    "    sum_density = 0\n",
    "    \n",
    "    for item in xrange(nums_neigh):\n",
    "        \n",
    "        sum_density += lrd_list[item]/lrd_list[index]\n",
    "        \n",
    "        #sum_density += lrd_list[item]\n",
    "    return sum_density/nums_neigh\n",
    "    #return sum_density/(nums_neigh*lrd_list[index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# 计算第k可达距离\n",
    "k_distance, neighbor_info = k_Distance(data,3) \n",
    "\n",
    "# 计算局部可达密度 local reachability density \n",
    "lrd_list = reach_density(data, k_distance, neighbor_info) \n",
    "\n",
    "# 样本点数量\n",
    "nums = data.shape[0] \n",
    "\n",
    "# 最终的LOF因子列表\n",
    "lof_list = []\n",
    "\n",
    "for index in xrange(nums):\n",
    "    \n",
    "    lof = cal_lof(index, data, neighbor_info, lrd_list)\n",
    "    lof_list.append(lof)\n",
    "\n",
    "# 设置阈值，大于阈值的认为是离群值\n",
    "threshold = 2\n",
    "boolean_array = [outlier_item > threshold for outlier_item in lof_list]\n",
    "\n",
    "indicy = []\n",
    "for key, value in enumerate(boolean_array):\n",
    "    if value:\n",
    "        indicy.append(key)\n",
    "print indicy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lof_list[73]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3 在二维图上表示每个样本点的LOF值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 人为产生离群值，检测算法\n",
    "#test_2d.loc[100]={'mileage':500,'price':1500}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = test_2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 设置阈值，大于阈值的认为是离群值\n",
    "threshold = 2\n",
    "boolean_array = [outlier_item > threshold for outlier_item in lof]\n",
    "\n",
    "indicy = []\n",
    "for key, value in enumerate(boolean_array):\n",
    "    if value:\n",
    "        indicy.append(key)\n",
    "print indicy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(7,7))\n",
    "\n",
    "for index,size in zip(xrange(nums),lof):\n",
    "    if index in indicy:\n",
    "        plt.scatter(data['mileage'][index], \n",
    "                    data['price'][index], \n",
    "                    s=np.exp(lof[index])*50,\n",
    "                    c='#efab40',\n",
    "                    alpha=0.6,\n",
    "                    marker='o')\n",
    "        plt.text(data['mileage'][index]-np.exp(lof[index])*50, \n",
    "                 data['price'][index]-np.exp(lof[index])*50,\n",
    "                str(round(lof[index],2)))\n",
    "    else:\n",
    "        plt.scatter(data['mileage'][index], \n",
    "                    data['price'][index], \n",
    "                    s=np.exp(lof[index])*50,\n",
    "                    c='#5dbe80',\n",
    "                    alpha=0.6,\n",
    "                    marker='o')\n",
    "        plt.text(data['mileage'][index]-np.exp(lof[index])*50, \n",
    "                 data['price'][index]-np.exp(lof[index])*50,\n",
    "                str(round(lof[index],2)),\n",
    "                fontsize=7)\n",
    "\n",
    "plt.xlabel('mileage')\n",
    "plt.ylabel('price')\n",
    "plt.grid('off')\n",
    "#sns.set_style('white')\n",
    "plt.savefig('./outlier.eps',format='eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 练习：尝试使用 oneclassSVM、IsolationForest和Robust covariance方法对离群价格进行检测？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><font color=\"#888\">Copyright &copy; 2017 Sid Ian@Boya, <a href=\"mailto:xiaodongy@boyabigdata.cn\">xiaodongy@boyabigdata.cn</a>. All Rights Reserved.</font></center>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "333px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
